2024 USRA Data Analysis

Introduction

The following tracks the data analysis for our 2024 Undergraduate Student Research Award (USRA) project, “What accounts for the relationship between disgust and religiosity?” This is the master analysis document, so all analyses are documented here—with any other code called within a chunk below. The goals of this analysis are to (1) clean the data, (2) check for internal reliability of our measures, (3) generate descriptive statistics, (4) visualize the distribution of our variables, (5) test our hypotheses, and (6) conduct exploratory analyses as appropriate.

Hypotheses

Individual differences in disgust sensitivity have repeatedly been shown to relate to religiosity, as predicted by the behavioral immune system model of disgust, leading to some scholars suggesting that religiosity is an evolved disease avoidance strategy. However, there are multiple mechanisms that aim to account for how religion contributes to disease avoidance, and these different mechanisms provide predictions for the psychological mediators of the relationship between disgust and religiosity. This study aims to test three such mediators, to better differentiate and understand why disgust sensitivity relates to religiosity. These mediators are represented by the following hypotheses.

1. The relationship between disgust sensitivity and religiosity is motivated by adherence to traditional practices.

   • Strong version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be fully mediated by conventionalism.

   • Weak version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be partially mediated by conventionalism.
2. The relationship between disgust sensitivity and religiosity is motivated by out-group avoidance.

   • Strong version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be fully mediated by ethnocentrism.

   • Weak version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be partially mediated by ethnocentrism.
3. The relationship between disgust sensitivity and religiosity is motivated by a monogamous mating strategy.

   • Strong version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be fully mediated by restricted sociosexual attitudes.

   • Weak version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be partially mediated by restricted sociosexual attitudes.

We will test these hypotheses through mediation analysis, the specifics of which will be reviewed in the Testing Hypotheses section.

Data Cleaning

All data files are saved in .csv format in the data directory, and the steps that were taken to process these files are documented in data/data_preparation_notes.txt. To clean the data for analysis, we will call a pre-written script (code/data_cleaning.R) that returns a data frame called data. The script also saves this data frame as data/data_clean.csv.

# Call the data cleaning script
source("./code/data_cleaning.R")

Warning: package 'groundhog' was built under R version 4.4.1

Attached: 'Groundhog' (Version: 3.2.1)

Tips and troubleshooting: https://groundhogR.com

Note: The Three Domains of Disgust Scale (TDDS) is generally administered with seven response options (Olatunji et al., 2012), only anchored with 0 (Not disgusting at all) and 6 (Extremely disgusting). Due to a survey coding error, our sample completed the TDDS with only six response options and the same anchors. Because of this, each TDDS item has a range from 0 to 5 (rather than 6), and raw scores for the full scale and subscales will be underestimated accordingly.

Variables in the Analysis

The data data frame contains the following variables (and column names).

• Participant ID (ID):

  • A unique randomly generated 10 digit string of numbers that identifies each participant.

• Start Time (start_time) and Finish Time (finish_time):

  • The date and time (down to the minute in PST) the participant started and finished the survey in POSIXct format (YYYY-MM-DD HH:MM:SS).

• Age (age):

  • The age of the participant in years.

• Sex (sex):

  • The sex of the participant as a factor variable (“Male” or “Female”).

• Gender (gender):

  • The gender of the participant. Responses were standardized into six categories formulated based on the genders reported. To see which original responses were put into these bins, see the appropriate section of the code/data_cleaning.R script. The categories are as follows:

    1. Male
    2. Female
    3. Non-binary
    4. Gender Queer
    5. Gender Fluid
    6. Agender

• Ethnicity (ethnicity):

  • The ethnicity of the participant. Participants were provided with a list of potential ethnicities to choose from as well as the options of writing in their own or “Rather Not Say”. Those who reported multiple ethnicities were categorized as “Mixed”. There are 11 unique ethnicities for this sample:

    1. South Asian
    2. East Asian
    3. South East Asian
    4. Latino or Hispanic
    5. White or Caucasian
    6. White or Sapharic Jew
    7. African
    8. Black or African American
    9. Black British
    10. Middle Eastern
    11. Mixed

• Nationality (nationality):

  • The nationality of the participant. Participants were provided with a list of potential ethnicities to choose from as well as the option of writing in their own. There were 41 unique nationalities reported by participants.

• Religious Affiliation (religious_affiliation):

  • The religion of participants. Participants were provided with a list of potential religions to choose from as well as the option of writing in their own. There were 12 religions reported, including “None”. Some participants wrote in “Catholic” or “Roman Catholic”, and were categorized as Christian (their Christian Affiliation was imputed as “Catholic”).

• Christian Affiliation (christian_affiliation):

  • The Christian sect or denomination reported by participants. Participants were provided with a list of potential groups to choose from as well as the option of writing in their own. Including “None” and “Non-religious Christian”, there were 16 affiliations reported.

• Religious Importance (religious_importance):

  • Participants were instructed to “Use the slider to indicate how important religion is in your life from 0 (not important at all) to 100 (extremely important).” This single-item measure will be used as a secondary indicator of religiosity, to examine whether results generalize beyond our primary DV of the Centrality of Religiosity Scale.

• The Centrality of Religiosity Scale (CRS):

  • The Centrality of Religiosity Scale (CRS-5) is a short measure of religiosity intended to triangulate the measurement of the multidimensional concept of religiosity by getting at five dimensions: intellect, ideology, public practice, private practice, and experience (Huber & Huber, 2012). Recent updates ensure that this measure is suitable for a cross-cultural sample. This measure shows strong concurrent validity with other religiosity measures and strong internal reliability.
  • Because the responses to items of the CRS-5 are qualitatively distinct, there are different response options, but the response format is the same, with a range of options that are scored from 1-5 (see the specific response option below their respective items). Where it is possible to get objective quantities for frequency (i.e., for religious service and prayer), the response options are the same. There are seven response options for these items, and they are collapsed into five scoring options (see below). Where objective quantities are potentially less frequent and are more difficult to capture (i.e., thinking about religious issues), the response options range from never to very often. Where objective quantities are not possible (i.e., degree of belief in God or the divine), response options are captured on a scale from not at all to very much so. o Once responses are scored according to the 1-5 system, all items are summed and divided by the number of items (5) to provide a total score ranging from 1.0-5.0.

• The Three Domains of Disgust Scale (full, TDDS_f; pathogen, TDDS_p; sexual, TDDS_s, moral, TDDS_m):

  • The TDDS was developed by Tybur et al. (2009) in response to a lack of theoretical and empirical grounding for prior measures of disgust sensitivity. Based on evolutionary theory, they predicted that items which people find disgusting should form three factors, one for pathogen, sexual, and moral disgust. Across multiple studies, they demonstrate that a three-factor solution is most parsimonious and fits the data well, that the full scale has good concurrent validity, and that the subscales have good convergent and divergent validity. Subsequent investigations have confirmed the three factor structure of the scale, although they find poorer validity evidence for the moral disgust subscale (Olatunji et al., 2012). All three subscales will be used in exploratory analyses and bivariate correlations, but the pathogen disgust scale will be used to test our central three hypotheses and six specific statistical hypotheses.

  • Scores for the full scale and each subscale were calculated by summing item values.

• The Revised Sociosexual Orientation Inventory (full, SOI_f; behavior, SOI_b; attitudes, SOI_a; desire, SOI_d):

  • The original Sociosexual Orientation Inventory (SOI) was a 7-item global measure of sociosexual orientation (Simpson & Gangestad, 1991), but this measure included behavioral, attitudinal, and (one) desire items together, with different response formats. This is problematic because each of these three components are differentiable both conceptually and empirically, and the different response formats often lead to an improper weighting of each of the three components contributing to the total score, resulting in poor reliability. Penke & Asendorpf (2008) developed the SOI-R to address these, and other, problems. The SOI-R is a 9-item measure with three items each tailored towards behaviors, attitudes, and desires, which they differentiated based on exploratory and confirmatory factor analysis as well as convergent and divergent validation. Each of the subscales show adequate reliability, along with the global score. As we are interested in operationalizing a monogamous mating strategy, it is participants intentions that we were most interested in measuring. In this case, one’s desire and behavior may, or may not, correspond to their intentions, which would best be operationalized as sociosexual attitudes. Therefore, we will use the attitudes subscale of the SOI-R to operationalize a monogamous mating strategy. This will be used to test our mediational hypothesis involving monogamous sexual strategies.

  • Scores for the full scale and each subscale were calculated by taking the arithmetic mean of item values.

• The In-Group Preference Subscale of the Short Form of the Generalized Ethnocentrism Scale (SFGENE-7; full, GENE_f; preference, GENE_p; superiority, GENE_s):

  • The SFGENE-7 is simply a short form of the Generalized Ethnocentrism Scale (Neuliep, 2002; Neuliep & McCroskey, 1997). The GENE is composed of 22 items, 15 of which are scored, and out of all the measures that we have, the GENE would be by far the longest. By reducing the number of items to 7 with the SFGENE-7, we can allow for the recruitment of more participants. The SFGENE-7 has shown adequate reliability despite a loss of items, as well as evidence for convergent validity with constructs like tolerance and multicultural ideology. In addition, it has been broken down via factor analysis into two distinct subscales, the first three items comprising the in-group preference scale and the last four items comprosing the in-group superiority scale. Each of these subscales also have strong internal reliabilities, and the in-group preference scale shows stronger negative correlations with tolerance and multicultural identity than the in-group superiority scale. We will use use only the first three items to measure the in-group preference aspect of ethnocentrism in order to operationalize out-group avoidance, particularly because the items comprising the in-group superiority subscale are quite similar at face value to items in the conventionalism scale that we are using to measure the tendency towards adherence to tradition. The in-group preference subscale will be used to test the mediational hypothesis involving out-group avoidance.

  • Items were summed to provide the total, preference, and superiority scores.

• Conventionalism Subscale of Aggression-Submission-Conventionalism Scale of Authoritarianism (ASC-C; CONV) and additional Traditionalism Items (TRAD; together, CONV_f):

  • Traditionalism is often operationalized as the traditionalism subscale of the Right Wing Authoritarianism (RWA) scale (Duckitt et al., 2010), but there are problems with this operationalization for our purposes. The RWA scale, and the traditionalism scale in particular, uses explicitly religious, sexual, and, some have argued, prejudicial content. Dunwoody & Funke (2016) address these limitations of the RWA scale by developing the Aggression-Submission-Conventionalism scale of Right Wing Authoritarianism that omits any politically charged, evaluative, prejudicial, or explicitly religious or sexual language. This removes the problems of tautology when predicting, in our case, religiosity, while retaining the same factor structure and similar reliability and validity as previous RWA scales. The conventionalism subscale, in particular, is meant to measure the same thing as the traditionalism subscale of the RWA scale but without religious language. Despite removing religious language, it retains a positive relationship with religiosity, and it seems to have better face validity than the traditionalism subscale for measuring adherence to traditional norms/practices. It is also quite short, allowing for quick administration without sacrificing psychometric properties.
  • In an attempt to improve our measure of traditionalism, we added some items which hopefully will get at our construct more specifically. These items are the TRAD items in the data frame. We will calculate Chronbach’s alpha and only accept these items in the measure if we keep an alpha above .8. If the alpha is below .8, we will remove items until it performs at this level.
  • Items were summed to provide each of the total scores.

Internal Reliability

To quickly assess the internal reliability of our scales, we will calculate Chronbach’s alpha for each of our scales. We will start with the CONV_f scale, because we may remove items we added to it.

# Defining the packages to use to do the alpha analysis
pkg <- c("tidyverse", "psych")
  # tidyverse if not attached from data cleaning script (for pipe and select function)
  # psych for alpha (and later descriptives)

# Loading the groundhog package to install packages from a certain date
library(groundhog)

# Reading in the packages with the groundhog package for reproducibility
  # Message suppressed in order to allow for rendering of quarto document
suppressMessages(groundhog.library(pkg = pkg, date = "2024-06-01"))

# Remove the pkg vector
rm(pkg)

Traditionalism

# Chronbach's Alpha for CONV_f (which includes all CONV and TRAD items)
alpha(x = data %>% select(starts_with("CONV."), starts_with("TRAD.")))

Reliability analysis   
Call: alpha(x = data %>% select(starts_with("CONV."), starts_with("TRAD.")))

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
      0.83      0.83    0.87      0.29 4.9 0.015    3 0.86     0.28

    95% confidence boundaries 
         lower alpha upper
Feldt      0.8  0.83  0.86
Duhachek   0.8  0.83  0.86

 Reliability if an item is dropped:
       raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CONV.1      0.82      0.82    0.86      0.30 4.6    0.015 0.034  0.29
CONV.2      0.80      0.80    0.84      0.27 4.0    0.017 0.031  0.25
CONV.3      0.82      0.82    0.85      0.29 4.6    0.015 0.029  0.30
CONV.4      0.81      0.81    0.85      0.28 4.4    0.016 0.034  0.28
CONV.5      0.82      0.81    0.85      0.28 4.4    0.016 0.034  0.28
CONV.6      0.82      0.82    0.85      0.29 4.5    0.016 0.028  0.29
TRAD.1      0.84      0.84    0.87      0.33 5.4    0.014 0.023  0.31
TRAD.2      0.82      0.82    0.85      0.29 4.5    0.016 0.030  0.28
TRAD.3      0.82      0.82    0.86      0.29 4.5    0.015 0.034  0.28
TRAD.4      0.82      0.81    0.85      0.28 4.4    0.016 0.031  0.29
TRAD.5      0.81      0.80    0.84      0.27 4.1    0.017 0.028  0.26
TRAD.6      0.81      0.81    0.85      0.28 4.3    0.016 0.030  0.28

 Item statistics 
         n raw.r std.r r.cor r.drop mean  sd
CONV.1 289  0.53  0.53  0.47   0.42  2.3 1.4
CONV.2 289  0.76  0.75  0.74   0.68  3.4 1.5
CONV.3 289  0.54  0.55  0.52   0.43  3.4 1.5
CONV.4 289  0.64  0.63  0.58   0.53  2.7 1.5
CONV.5 289  0.63  0.63  0.59   0.54  2.2 1.5
CONV.6 289  0.57  0.58  0.55   0.47  3.9 1.3
TRAD.1 289  0.27  0.27  0.17   0.14  1.5 1.4
TRAD.2 289  0.57  0.57  0.53   0.47  3.9 1.4
TRAD.3 289  0.57  0.57  0.51   0.46  2.0 1.5
TRAD.4 289  0.63  0.63  0.60   0.53  3.0 1.5
TRAD.5 289  0.73  0.73  0.73   0.66  3.7 1.4
TRAD.6 289  0.63  0.64  0.61   0.54  3.7 1.4

Non missing response frequency for each item
          0    1    2    3    4    5    6 miss
CONV.1 0.13 0.17 0.27 0.25 0.12 0.03 0.02    0
CONV.2 0.04 0.07 0.14 0.28 0.24 0.13 0.10    0
CONV.3 0.04 0.06 0.12 0.29 0.25 0.15 0.09    0
CONV.4 0.07 0.13 0.27 0.26 0.13 0.08 0.06    0
CONV.5 0.13 0.18 0.31 0.21 0.10 0.04 0.03    0
CONV.6 0.01 0.03 0.09 0.25 0.30 0.19 0.12    0
TRAD.1 0.31 0.25 0.22 0.15 0.03 0.02 0.01    0
TRAD.2 0.02 0.04 0.10 0.20 0.31 0.17 0.15    0
TRAD.3 0.19 0.23 0.25 0.16 0.11 0.04 0.01    0
TRAD.4 0.06 0.09 0.23 0.25 0.23 0.09 0.06    0
TRAD.5 0.03 0.05 0.09 0.26 0.32 0.13 0.12    0
TRAD.6 0.02 0.07 0.08 0.25 0.30 0.20 0.08    0

• The raw alpha value (.831; standardized .830) is acceptable here, so we will use our full CONV_f variable to operationalize traditionalism.

• To compare these values to those for the original conventionalism measure, I will calculate its’ alpha as well.

# Chronbach's Alpha for CONV (only CONV items)
alpha(x = data %>% select(starts_with("CONV.")))

Reliability analysis   
Call: alpha(x = data %>% select(starts_with("CONV.")))

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
      0.74      0.74    0.76      0.33 2.9 0.024    3 0.97     0.38

    95% confidence boundaries 
         lower alpha upper
Feldt      0.7  0.74  0.79
Duhachek   0.7  0.74  0.79

 Reliability if an item is dropped:
       raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CONV.1      0.73      0.73    0.75      0.36 2.8    0.025 0.028  0.40
CONV.2      0.67      0.67    0.70      0.29 2.1    0.031 0.038  0.21
CONV.3      0.72      0.72    0.70      0.34 2.6    0.025 0.018  0.39
CONV.4      0.70      0.70    0.71      0.32 2.4    0.028 0.031  0.36
CONV.5      0.69      0.69    0.70      0.31 2.3    0.029 0.030  0.35
CONV.6      0.72      0.71    0.70      0.33 2.5    0.026 0.020  0.39

 Item statistics 
         n raw.r std.r r.cor r.drop mean  sd
CONV.1 289  0.59  0.59  0.45   0.39  2.3 1.4
CONV.2 289  0.75  0.75  0.67   0.59  3.4 1.5
CONV.3 289  0.62  0.63  0.56   0.43  3.4 1.5
CONV.4 289  0.68  0.67  0.58   0.50  2.7 1.5
CONV.5 289  0.70  0.70  0.62   0.53  2.2 1.5
CONV.6 289  0.63  0.64  0.58   0.45  3.9 1.3

Non missing response frequency for each item
          0    1    2    3    4    5    6 miss
CONV.1 0.13 0.17 0.27 0.25 0.12 0.03 0.02    0
CONV.2 0.04 0.07 0.14 0.28 0.24 0.13 0.10    0
CONV.3 0.04 0.06 0.12 0.29 0.25 0.15 0.09    0
CONV.4 0.07 0.13 0.27 0.26 0.13 0.08 0.06    0
CONV.5 0.13 0.18 0.31 0.21 0.10 0.04 0.03    0
CONV.6 0.01 0.03 0.09 0.25 0.30 0.19 0.12    0

• As you can see, the CONV_f variable has stronger internal reliability than the CONV items alone.

Religiosity

# Chronbach's Alpha for Centrality of Religiosity Items
alpha(x = data %>% select(starts_with("CRS.")))

Reliability analysis   
Call: alpha(x = data %>% select(starts_with("CRS.")))

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean  sd median_r
      0.92      0.92    0.91      0.69  11 0.0069  2.6 1.3     0.72

    95% confidence boundaries 
         lower alpha upper
Feldt     0.90  0.92  0.93
Duhachek  0.91  0.92  0.93

 Reliability if an item is dropped:
      raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
CRS.1      0.93      0.93    0.92      0.78 13.8   0.0064 0.002  0.80
CRS.2      0.89      0.89    0.88      0.68  8.4   0.0094 0.012  0.66
CRS.3      0.90      0.90    0.88      0.69  8.8   0.0090 0.016  0.70
CRS.4      0.89      0.89    0.86      0.66  7.7   0.0104 0.011  0.65
CRS.5      0.89      0.89    0.87      0.67  8.1   0.0096 0.015  0.65

 Item statistics 
        n raw.r std.r r.cor r.drop mean  sd
CRS.1 289  0.73  0.76  0.65   0.63  2.6 1.1
CRS.2 289  0.90  0.89  0.87   0.83  3.0 1.5
CRS.3 289  0.88  0.88  0.84   0.81  2.3 1.4
CRS.4 289  0.93  0.92  0.91   0.87  2.5 1.7
CRS.5 289  0.90  0.90  0.88   0.85  2.5 1.4

Non missing response frequency for each item
         1    2    3    4    5 miss
CRS.1 0.14 0.37 0.28 0.13 0.08    0
CRS.2 0.18 0.29 0.15 0.11 0.27    0
CRS.3 0.42 0.22 0.13 0.09 0.14    0
CRS.4 0.46 0.15 0.08 0.03 0.28    0
CRS.5 0.33 0.26 0.16 0.11 0.14    0

• The CRS shows strong internal consistency here (alpha = .919; standardized = .919).

Pathogen Disgust

# Chronbach's Alpha for Pathogen Disgust Items
alpha(x = data %>% select(., c("TDDS.12", "TDDS.15", "TDDS.9", "TDDS.3", "TDDS.21", "TDDS.18", "TDDS.6")))

Reliability analysis   
Call: alpha(x = data %>% select(., c("TDDS.12", "TDDS.15", "TDDS.9", 
    "TDDS.3", "TDDS.21", "TDDS.18", "TDDS.6")))

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
      0.77      0.77    0.76      0.33 3.4 0.021    3 0.91     0.34

    95% confidence boundaries 
         lower alpha upper
Feldt     0.73  0.77  0.81
Duhachek  0.73  0.77  0.81

 Reliability if an item is dropped:
        raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
TDDS.12      0.74      0.74    0.72      0.33 2.9    0.023 0.0094  0.34
TDDS.15      0.72      0.72    0.70      0.30 2.6    0.026 0.0063  0.28
TDDS.9       0.72      0.73    0.70      0.31 2.7    0.025 0.0064  0.31
TDDS.3       0.75      0.75    0.74      0.34 3.1    0.023 0.0091  0.34
TDDS.21      0.75      0.75    0.74      0.33 3.0    0.023 0.0114  0.35
TDDS.18      0.75      0.75    0.74      0.34 3.1    0.022 0.0095  0.35
TDDS.6       0.75      0.76    0.74      0.34 3.1    0.023 0.0076  0.34

 Item statistics 
          n raw.r std.r r.cor r.drop mean  sd
TDDS.12 289  0.65  0.65  0.57   0.49  3.4 1.4
TDDS.15 289  0.73  0.74  0.70   0.61  3.3 1.3
TDDS.9  289  0.71  0.71  0.67   0.57  2.5 1.4
TDDS.3  289  0.58  0.61  0.51   0.44  4.0 1.1
TDDS.21 289  0.65  0.63  0.52   0.47  3.0 1.5
TDDS.18 289  0.63  0.61  0.51   0.44  3.1 1.5
TDDS.6  289  0.61  0.60  0.50   0.44  1.9 1.4

Non missing response frequency for each item
           0    1    2    3    4    5 miss
TDDS.12 0.02 0.09 0.15 0.21 0.25 0.28    0
TDDS.15 0.01 0.11 0.15 0.26 0.24 0.24    0
TDDS.9  0.07 0.20 0.28 0.21 0.15 0.10    0
TDDS.3  0.00 0.03 0.09 0.17 0.23 0.48    0
TDDS.21 0.06 0.15 0.18 0.22 0.15 0.24    0
TDDS.18 0.04 0.14 0.18 0.19 0.17 0.27    0
TDDS.6  0.17 0.24 0.26 0.18 0.08 0.06    0

• The TDDS Pathogen Disgust (TDDS_p) shows adequate internal reliability here (alpha = .769; standardized = .772).

Sociosexual Attitudes

# Chronbach's Alpha for the Attitudes Subscale of the SOI-R
alpha(x = data %>% select(., c("SOI.4", "SOI.5", "SOI.6")))

Reliability analysis   
Call: alpha(x = data %>% select(., c("SOI.4", "SOI.5", "SOI.6")))

  raw_alpha std.alpha G6(smc) average_r S/N  ase mean  sd median_r
      0.81      0.81    0.74      0.58 4.2 0.02  5.2 2.4      0.6

    95% confidence boundaries 
         lower alpha upper
Feldt     0.76  0.81  0.84
Duhachek  0.77  0.81  0.84

 Reliability if an item is dropped:
      raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
SOI.4      0.76      0.76    0.62      0.62 3.2    0.028    NA  0.62
SOI.5      0.69      0.69    0.52      0.52 2.2    0.037    NA  0.52
SOI.6      0.75      0.75    0.60      0.60 3.0    0.029    NA  0.60

 Item statistics 
        n raw.r std.r r.cor r.drop mean  sd
SOI.4 289  0.83  0.83  0.70   0.63  5.9 2.7
SOI.5 289  0.88  0.87  0.78   0.70  4.3 2.9
SOI.6 289  0.84  0.84  0.71   0.64  5.3 2.8

Non missing response frequency for each item
         1    2    3    4    5    6    7    8    9 miss
SOI.4 0.11 0.04 0.08 0.06 0.11 0.09 0.15 0.11 0.25    0
SOI.5 0.29 0.11 0.10 0.04 0.10 0.09 0.09 0.06 0.13    0
SOI.6 0.17 0.04 0.11 0.05 0.14 0.09 0.08 0.12 0.19    0

• Despite only comprising three items, the attitudes subscale of the SOI-R (SOI_a) seems to have good internal consistency here (alpha = .806; standardized = .806).

Preferences Subscale of SFGENE-7

# Chronbach's Alpha for the Preferences Subscale of the GENE
alpha(x = data %>% select(., c("GENE.1", "GENE.2", "GENE.3")))

Reliability analysis   
Call: alpha(x = data %>% select(., c("GENE.1", "GENE.2", "GENE.3")))

  raw_alpha std.alpha G6(smc) average_r S/N  ase mean   sd median_r
      0.52      0.53    0.55      0.27 1.1 0.05  2.3 0.68    0.095

    95% confidence boundaries 
         lower alpha upper
Feldt     0.42  0.52  0.61
Duhachek  0.43  0.52  0.62

 Reliability if an item is dropped:
       raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
GENE.1     0.801     0.804   0.672     0.672 4.09    0.023    NA 0.672
GENE.2     0.093     0.093   0.049     0.049 0.10    0.107    NA 0.049
GENE.3     0.174     0.174   0.095     0.095 0.21    0.097    NA 0.095

 Item statistics 
         n raw.r std.r r.cor r.drop mean   sd
GENE.1 289  0.54  0.53 0.096  0.078  3.2 0.95
GENE.2 289  0.81  0.82 0.766  0.536  1.7 0.89
GENE.3 289  0.80  0.80 0.734  0.472  1.9 0.99

Non missing response frequency for each item
          1    2    3    4    5 miss
GENE.1 0.05 0.12 0.47 0.27 0.09    0
GENE.2 0.51 0.32 0.12 0.03 0.01    0
GENE.3 0.47 0.29 0.18 0.06 0.01    0

• The internal consistency of the preferences subscale of the SFGENE-7 is quite low (alpha = .523; standardized = .528). However, the removal of Item 1 improves the alpha statistic drastically (to alpha = .801; standardized = .804). Because of this, we will run the analysis with both the full Preferences Subscale and only Items 2 and 3. I will calculate the scores for the latter variable and call it GENE_p_2.

# Creating the GENE_p_2 variable as the sum of items 2 and 3
data <- data %>% 
  mutate(
    GENE_p_2 = rowSums(select(., c(GENE.2, GENE.3)))
  )

Descriptive Statistics

Sample Size and Time Taken

First, we will see how many participants we have and how long it took participants to complete the survey.

# Calculate the number of participants in the sample
nrow(data)

[1] 289

# Calculate descriptive statistics for the difference between the finish and start time (in minutes)
describe(as.numeric(difftime(data$finish_time, data$start_time, units = "mins")))

   vars   n  mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 289 14.39 58.1      8    8.94 2.97   3 821   818 12.25   154.08 3.42

• In total, we have N = 289 participants with complete data.

• Although the mean time taken is around 14.39, the SD is huge, and the maximum is a massive 821 minutes (over 13 hours). The trimmed mean of 8.94 is therefore a better estimate.

Categorical Variables

Next, we will describe the sample in terms of the categorical variables. The following generates frequency tables for all categorical variables.

# Create a frequency table for each factor variable
lapply(data[sapply(data, is.factor)], table)

$sex

  Male Female 
   146    143 

$gender

     Agender       Female Gender Fluid Gender Queer         Male   Non-binary 
           2          135            1            1          144            4 

$ethnicity

                          African         Black or African American 
                               52                                 4 
                        Caribbean                        East Asian 
                                0                                14 
               Latino or Hispanic                    Middle Eastern 
                               62                                 4 
                            Mixed Native American or Alaskan Native 
                               11                                 0 
                      South Asian                White or Causasian 
                                8                               122 
            White or Sapharic Jew                     Black British 
                                6                                 1 
                    White Mexican               Romani or Traveller 
                                0                                 0 
                 South East Asian                    Rather Not Say 
                                2                                 3 

$nationality

       Algeria      Australia        Austria        Belgium         Brazil 
             1              7              1              1              3 
        Canada          Chile          China       Columbia Czech Republic 
            30             13              1              2              5 
         Egypt        Eritrea        Finland         France        Germany 
             1              1              2              3              4 
         Ghana         Greece        Hungary          India        Ireland 
             1              5              3              1              1 
        Israel          Italy          Japan          Kenya         Mexico 
             2              5              1              6             39 
   New Zealand         Poland       Portugal        Romania   Saudi Arabia 
             6             14             25              1              1 
      Slovakia   South Africa          Spain         Sweden        Tunisia 
             1             46              7              5              1 
        Turkey United Kingdom  United States      Venezuela        Vietnam 
             2             21             17              1              1 
      Zimbabwe 
             1 

$religious_affiliation

      Agnostic Anti-religious        Atheist       Buddhist      Christian 
             6              1              2              2            123 
         Deist          Hindu         Muslim           None          Pagan 
             1              1             11            139              1 
     Spiritual   Spiritualist 
             1              1 

$christian_affiliation

               Anglican               Apostolic                 Baptist 
                      1                       1                      11 
               Catholic       Congregationalist             Evangelical 
                     68                       1                       4 
      Jehovah’s Witness                Lutheran               Methodist 
                      1                       6                       9 
Non-religious Christian                    None                Orthodox 
                      1                       4                       4 
            Pentecostal            Presbyterian              Protestant 
                      1                       2                       8 
  Seventh Day Adventist 
                      1 

The output is likely difficult to read, so I will break it down here.

• Sex:

  • Due to our use of the balance sex tool on Prolific, we have a very similar number of males (n = 146, 50.5%) and females (n = 143, 49.5%).

• Gender:

  • Male: 144 (49.8%)

  • Female: 135 (46.7%)

  • Non-Binary: 4 (1.4%)

  • Agender: 2 (.7%)

  • Gender Fluid: 1 (.3%)

  • Gender Queer: 1 (.3%)

• Ethnicity:

  • The sample is predominantly White or Caucasian (n = 122, 42.2%), followed by Latino or Hispanic (n = 62, 21.5%), African (n = 52, 18.0%), East Asian (n = 14, 4.8%), and Mixed (n = 11, 3.8%). All other ethnic categories have n = 8 or fewer participants in them (≤2.8%).

• Nationality:

  • There are a large number of nationalities represented in our sample. The most common nationality is South Africa (n = 46, 15.9%), followed by Mexico (n = 39, 13.5%), Canada (n = 30, 10.4%), Portugal (n = 25, 8.7%), the United Kingdom (n = 21, 7.3%), and the United States (n = 17, 5.9%). All other nationalities have 14 or fewer participants (≤4.8%).

• Religious Affiliation:

  • The modal response for Religious Affiliation in this sample is “None” (n = 139, 48.1%). The most common religious affiliation other than this is, by far, Christian (n = 123, 42.6%). The next most common categories is Muslim (n = 11, 3.8%) and Agnostic (n = 6, 2.1%). All other categories have two or fewer participants (≤.7%).

• Christian Affiliation:

  • Of those who identified as Christian, most are Catholic (n = 68, 55.3%). There are also a substantial number of Baptists (n = 11, 8.9%), Methodists (n = 9, 7.3%), and Protestants (n = 8, 6.5%). All other categories have four or fewer participants (≤3.3%).

Numeric Variables

Now we will describe our numeric variables starting with our only demographic numeric variable, age.

# Generate descriptive statistics for age
describe(data$age)

   vars   n  mean    sd median trimmed  mad min max range skew kurtosis   se
X1    1 289 30.85 10.85     28   29.18 7.41  18  86    68 1.74     3.84 0.64

• The average age in our sample is 30.85 (SD = 10.85), with a trimmed mean of 29.18, indicating little skew. This indicates that our sample tends to be younger adults, but with quite a bit of variability.

Now we will describe the important scales for our analysis.

# Describe important scales in our analysis
describe(data %>% select("religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))

                     vars   n  mean    sd median trimmed   mad min max range
religious_importance    1 289 30.70 35.63  12.00   26.08 17.79   0 100   100
CRS                     2 289  2.59  1.26   2.20    2.50  1.19   1   5     4
TDDS_f                  3 289 60.12 16.05  59.00   59.82 16.31  12 100    88
TDDS_p                  4 289 21.26  6.35  22.00   21.30  5.93   4  35    31
SOI_f                   5 289  3.58  1.62   3.22    3.48  1.81   1   8     7
SOI_a                   6 289  5.16  2.40   5.00    5.19  2.47   1   9     8
GENE_p                  7 289  6.81  2.03   7.00    6.72  1.48   3  15    12
GENE_p_2                8 289  3.57  1.71   3.00    3.34  1.48   2  10     8
CONV_f                  9 289 35.55 10.30  36.00   35.78 10.38   0  65    65
                      skew kurtosis   se
religious_importance  0.87    -0.83 2.10
CRS                   0.57    -1.06 0.07
TDDS_f                0.17    -0.29 0.94
TDDS_p               -0.07    -0.49 0.37
SOI_f                 0.53    -0.59 0.10
SOI_a                -0.03    -1.07 0.14
GENE_p                0.50     0.36 0.12
GENE_p_2              1.02     0.64 0.10
CONV_f               -0.24     0.44 0.61

• Other than a bit of skew for religious_importance and the GENE_p_2, there does not seem to be a problem here. Variables will be inspected more closely when testing assumptions for statistical tests.

Before moving on, we will calculate correlations between each of these variables.

# Calculate correlations for numerical variables
corr.test(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "TDDS_s", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))

Call:corr.test(x = data %>% select("age", "religious_importance", 
    "CRS", "TDDS_f", "TDDS_p", "TDDS_s", "SOI_f", "SOI_a", "GENE_p", 
    "GENE_p_2", "CONV_f"))
Correlation matrix 
                       age religious_importance   CRS TDDS_f TDDS_p TDDS_s
age                   1.00                -0.07 -0.09   0.02  -0.08  -0.09
religious_importance -0.07                 1.00  0.92   0.46   0.22   0.52
CRS                  -0.09                 0.92  1.00   0.47   0.24   0.53
TDDS_f                0.02                 0.46  0.47   1.00   0.76   0.82
TDDS_p               -0.08                 0.22  0.24   0.76   1.00   0.47
TDDS_s               -0.09                 0.52  0.53   0.82   0.47   1.00
SOI_f                 0.12                -0.28 -0.27  -0.38  -0.16  -0.54
SOI_a                 0.14                -0.39 -0.41  -0.45  -0.22  -0.59
GENE_p               -0.06                 0.07  0.07   0.02   0.01   0.07
GENE_p_2              0.03                 0.02  0.00  -0.02  -0.03   0.04
CONV_f                0.08                 0.33  0.31   0.12   0.01   0.13
                     SOI_f SOI_a GENE_p GENE_p_2 CONV_f
age                   0.12  0.14  -0.06     0.03   0.08
religious_importance -0.28 -0.39   0.07     0.02   0.33
CRS                  -0.27 -0.41   0.07     0.00   0.31
TDDS_f               -0.38 -0.45   0.02    -0.02   0.12
TDDS_p               -0.16 -0.22   0.01    -0.03   0.01
TDDS_s               -0.54 -0.59   0.07     0.04   0.13
SOI_f                 1.00  0.88  -0.11    -0.10  -0.10
SOI_a                 0.88  1.00  -0.15    -0.11  -0.15
GENE_p               -0.11 -0.15   1.00     0.88   0.08
GENE_p_2             -0.10 -0.11   0.88     1.00   0.03
CONV_f               -0.10 -0.15   0.08     0.03   1.00
Sample Size 
[1] 289
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
                      age religious_importance  CRS TDDS_f TDDS_p TDDS_s SOI_f
age                  0.00                 1.00 1.00   1.00   1.00   1.00  1.00
religious_importance 0.21                 0.00 0.00   0.00   0.01   0.00  0.00
CRS                  0.15                 0.00 0.00   0.00   0.00   0.00  0.00
TDDS_f               0.79                 0.00 0.00   0.00   0.00   0.00  0.00
TDDS_p               0.15                 0.00 0.00   0.00   0.00   0.00  0.21
TDDS_s               0.14                 0.00 0.00   0.00   0.00   0.00  0.00
SOI_f                0.05                 0.00 0.00   0.00   0.01   0.00  0.00
SOI_a                0.01                 0.00 0.00   0.00   0.00   0.00  0.00
GENE_p               0.35                 0.20 0.25   0.75   0.81   0.21  0.06
GENE_p_2             0.66                 0.79 1.00   0.72   0.67   0.55  0.09
CONV_f               0.17                 0.00 0.00   0.04   0.93   0.03  0.10
                     SOI_a GENE_p GENE_p_2 CONV_f
age                   0.42   1.00     1.00   1.00
religious_importance  0.00   1.00     1.00   0.00
CRS                   0.00   1.00     1.00   0.00
TDDS_f                0.00   1.00     1.00   1.00
TDDS_p                0.01   1.00     1.00   1.00
TDDS_s                0.00   1.00     1.00   0.74
SOI_f                 0.00   1.00     1.00   1.00
SOI_a                 0.00   0.29     1.00   0.32
GENE_p                0.01   0.00     0.00   1.00
GENE_p_2              0.06   0.00     0.00   1.00
CONV_f                0.01   0.16     0.62   0.00

 To see confidence intervals of the correlations, print with the short=FALSE option

• For each correlation the sample size is N = 289.

• Noteworthy for this analysis, pathogen disgust is positively and significantly related to both religious_importance (r = .22, p < .01) and the CRS (r = .24, p < .01), but sexual disgust is more strongly related to both religious_importance (r = .52, p < .01) and the CRS (r = .53, p < .01). This replicates previous work and provides initial support for the sexual strategies hypothesis.

  • In addition, there is a very strong correlation between religious_importance and the CRS (r = .92, p < .01).

  • However, neither the GENE_p nor the GENE_p_2 are significantly related to either religious_importance (r = .07, p = .20; r = .02, p = .02; respectively) or the CRS (r = .07, p = .25; r = .00, p = .00; respectively), suggesting that it could not mediate the relationship between disgust sensitivity and religiosity.

  • Unlike the GENE variables, CONV_f is significantly related to both religious_importance (r = .33, p < .01) and the CRS (r = .31, p < .01).

• Next, the only potential mediator that is significantly related to pathogen disgust is the SOI_a:

  • SOI_a: r = -.22, p = .01

  • GENE_p: r = .01, p = 1

    • GENE_p_2: r = -.03, p = 1

  • CONV_f: r = .01, p = 1

  • This suggests that the only plausible mediator (here, that is) of the relationship between disgust and religiosity could be monogamous mating strategies.

Data Visualization

Now we will visualization the distribution of our variables by creating bar plots and histograms (for categorical and numerical variables, respectively).

Categorical Variables

# Loop through each categorical variable in the data frame and create a bar chart
for (var in names(data[, sapply(data, is.factor)])) {
  p <- ggplot(data[, sapply(data, is.factor)], aes(x = !!sym(var))) + 
    geom_bar() +
    theme_minimal() +
    labs(title = paste("Bar Chart of", var), x = var, y = "Count") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) # Adjust x axis labels to make them readable (not all over eachother)
  # Print each bar chart
  print(p)
  # Remove the plot object so it doesn't clutter the environment
  rm(p)
}

# Remove the var object
rm(var)

Numeric Variables

# Loop through numeric variable identified to produce a histogram
for (var in names(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))) {
  p <- ggplot(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"), aes(x = .data[[var]])) + 
    geom_histogram(bins = 10, fill = "steelblue", color = "black") +
    theme_minimal() +
    labs(title = paste("Histogram of", var), x = var, y = "Frequency")

  # Print each histogram
  print(p)
  # Remove the p object so it doesn't clutter up the environment
  rm(p)
}

# Remove the var object from the environment
rm(var)

• The distributions of some of these variables indicate that they are non-normally distributed. With such a large sample (N = 289), this should not be much of a problem, particularly because we will be using bootstrapped confidence intervals for much of the analysis.

Testing Hypotheses

As above, the following are our hypotheses.

1. The relationship between disgust sensitivity and religiosity is motivated by adherence to traditional practices.

   • Strong version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be fully mediated by conventionalism.

   • Weak version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be partially mediated by conventionalism.
2. The relationship between disgust sensitivity and religiosity is motivated by out-group avoidance.

   • Strong version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be fully mediated by ethnocentrism.

   • Weak version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be partially mediated by ethnocentrism.
3. The relationship between disgust sensitivity and religiosity is motivated by a monogamous mating strategy.

   • Strong version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be fully mediated by restricted sociosexual attitudes.

   • Weak version of statistical hypothesis:

     • The positive relationship between parasite disgust sensitivity and religiosity will be partially mediated by restricted sociosexual attitudes.

The three hypotheses are each broken down into two statistical hypotheses, which represent their strong and weak version. In the strong version there is full mediation, and in the weak version there is only partial mediation. We have defined and preregistered inference criteria for both of these versions of each hypothesis.

Strong Version Inference Criterion:

For each model (each hypothesis) we will use the “K” method outlined by Beribisky et al. (2020) to infer full mediation. That is, if the proportion of variance explained by the indirect effect compared to the total effect is 80% or higher, we will infer full mediation. This will provide support for the strong version of the respective hypothesis.

Weak Version Inference Criterion:

Where there is not full mediation, for each model (each hypothesis) if the 95% confidence interval for the indirect effect of parasite disgust sensitivity on religiosity through the mediator does not contain zero, then we will infer partial mediation. This will provide support for the weak version of the respective hypothesis.

Data Preparation

Before testing our hypotheses, we will standardize each variable. We will also reverse the scores of SOI_a before standardizing it, because higher scores currently represent more unrestricted sociosexual attitudes, whereas we would like higher scores to represent more restricted sociosexual attitudes (as a proxy for a monogamous mating strategy).

# Reversing the scoring for the SOI_a (saving it as SOI_a_r) and standardizing all variables
data <- data %>% 
  mutate(
    SOI_a_r = (max(data$SOI_a) + 1) - SOI_a,
    SOI_a_r_z = scale(SOI_a_r),
    CRS_z = scale(CRS),
    religious_importance_z = scale(religious_importance),
    GENE_p_z = scale(GENE_p),
    GENE_p_2_z = scale(GENE_p_2),
    CONV_f_z = scale(CONV_f),
    TDDS_p_z = scale(TDDS_p)
  )

# Ensure that the variables are plain numeric vectors
data$TDDS_p_z <- as.numeric(data$TDDS_p_z)
data$CONV_f_z <- as.numeric(data$CONV_f_z)
data$CRS_z <- as.numeric(data$CRS_z)
data$SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
data$religious_importance_z <- as.numeric(data$religious_importance_z)
data$GENE_p_z <- as.numeric(data$GENE_p_z)
data$GENE_p_2_z <- as.numeric(data$GENE_p_2_z)

# Because the boot function was throwing a fit about this earlier
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

We will also load the packages involved in the mediation analysis.

# Defining the packages to use to do the mediation analysis
pkg <- c("car", "mediation")
  # car for testing for multicollinearity and bootstrapping
  # mediation for testing the indirect effect with bootstrapped CIs and accessing

# Reading in the packages with the groundhog package for reproducibility
  # Message suppressed in order to allow for rendering of quarto document
suppressMessages(groundhog.library(pkg = pkg, date = "2024-06-01"))

# Dropping the pkg vector
rm(pkg)

Testing the Baseline Model

Before moving on, we will first test the baseline model where we regress the dependent variable (CRS) on the independent variable (TDDS_p), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions.

# Fitting the model
baseline_model <- lm(CRS_z ~ TDDS_p_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(baseline_model, which = 1)

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(baseline_model))

   

   # QQ-plot for normality of residuals
   qqnorm(residuals(baseline_model))
   qqline(residuals(baseline_model), col = "red")

   

   • The histogram of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just around -1 SD and around 2 SD. The two clusters of similar errors in prediction may be due to sex differences in disgust and religiosity. In exploratory analyses, we will add sex as a control variable.

   • Because our inferential criteria do not require parametric tests, we should not be in too much trouble with a non-normality of residuals. We were relying on a parametric t-test to assess whether pathogen disgust is related to religiosity in the first place, however. Because we have failed the assumption of normality of residuals, we will use a bootstrapped test to make our inference here as well.

3. Multicollinearity: Not applicable becuase only one predictor

Summarizing the Model

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(baseline_model, col = "black", lwd = 1)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_baseline_model <- Boot(baseline_model, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_baseline_model) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original  bootBias   bootSE   bootMed
(Intercept) -8.5501e-16 0.0015289 0.057618 0.0036089
TDDS_p_z     2.4299e-01 0.0010657 0.059478 0.2411963

confint(boot_baseline_model) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1197906 0.1067505
TDDS_p_z     0.1345033 0.3605490

# Removing the CRS_z and TDDS_p_z vectors in the environment that the Boot package required outside of the data.frame
rm(CRS_z, TDDS_p_z)

• We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and religiosity.

• Based on the estimates, the standardized coefficient for pathogen disgust is around .24. This aligns with estimates from a meta-analysis of this relationship (Yu et al., 2022).

  • In addition, this relationship is significant, as the BCIs do not contain zero (upper = .35, lower = .12).

This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.

Hypothesis 1: Adherence to Traditional Practices

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religiosity.

# Fitting the models
H1_M1 <- lm(CONV_f_z ~ TDDS_p_z, data = data)
H1_M2 <- lm(CRS_z ~ CONV_f_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_M1, which = 1) # For model 1

   

   plot(H1_M2, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the first model, but there is a trend towards more negative residual values at the higher fitted values of the model for model 2. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_M1)) # For model 1

   

   hist(residuals(H1_M2)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H1_M1))
   qqline(residuals(H1_M1), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H1_M2))
   qqline(residuals(H1_M2), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CONV_f_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Traditionalism",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M1, col = "black", lwd = 1)

# Summarizing the model
summary(H1_M1)

Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4517 -0.6429  0.0402  0.6287  2.8524 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.379e-17  5.893e-02   0.000    1.000
TDDS_p_z     5.260e-03  5.903e-02   0.089    0.929

Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared:  2.766e-05, Adjusted R-squared:  -0.003457 
F-statistic: 0.00794 on 1 and 287 DF,  p-value: 0.9291

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = .005, t(287) = .089, p = .929), just as with the correlation above.

  • Although this precludes the possibility that adherence to tradition mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$CONV_f_z, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Traditionalism",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M2, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2 <- Boot(H1_M2, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_M2) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original  bootBias  bootSE   bootMed
(Intercept) -8.5089e-16 0.0026934 0.05539 0.0029616
CONV_f_z     3.1443e-01 0.0036837 0.05224 0.3164595

confint(boot_H1_M2) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1128222 0.1044219
CONV_f_z     0.2036523 0.4178083

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z)

• We can see from the Scatter Plot that there is a positive relationship between traditionalism and religiosity.

• Based on the estimates, the standardized coefficient for traditionalism is around .31.

  • In addition, this relationship is significant, as the BCIs do not contain zero (upper = .42, lower = .20).

Full Mediation Model

The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H1_full <- lm(CRS_z ~ TDDS_p_z + CONV_f_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_full, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_full))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H1_full))
   qqline(residuals(H1_full), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H1_full)

   TDDS_p_z CONV_f_z 
   1.000028 1.000028 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full <- Boot(H1_full, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_full) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE  bootMed
(Intercept) -8.3817e-16  0.00170100 0.053887 0.001676
TDDS_p_z     2.4134e-01 -0.00064514 0.056527 0.240610
CONV_f_z     3.1316e-01  0.00297386 0.048923 0.316369

confint(boot_H1_full) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1106346 0.1026964
TDDS_p_z     0.1328538 0.3560156
CONV_f_z     0.2100775 0.4057880

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results <- mediate(model.m = H1_M1,
        model.y = H1_full,
        treat = "TDDS_p_z",
        mediator = "CONV_f_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H1_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.00165     -0.03393         0.04    0.96    
ADE             0.24134      0.12018         0.35  <2e-16 ***
Total Effect    0.24299      0.13035         0.35  <2e-16 ***
Prop. Mediated  0.00678     -0.17652         0.19    0.96    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religiosity.

  • The coefficients for pathogen disgust (.24) and traditionalism (.31) are virtually identical in this model to the models above where they are stand-alone predictors. This suggests that the unique effects of pathogen disgust and traditionalism on religiosity over-and-above each other are completely independent effects.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = .001, BCI = [-.035, .04]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism.

Hypothesis 2: Out-Group Avoidance

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religiosity.

# Fitting the models
H2_M1 <- lm(GENE_p_z ~ TDDS_p_z, data = data)
H2_M2 <- lm(CRS_z ~ GENE_p_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_M1, which = 1) # For model 1

   

   plot(H2_M2, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_M1)) # For model 1

   

   hist(residuals(H2_M2)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H2_M1))
   qqline(residuals(H2_M1), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H2_M2))
   qqline(residuals(H2_M2), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are relatively normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Out-Group Avoidance (GENE_p)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M1, col = "black", lwd = 1)

# Summarizing the model
summary(H2_M1)

Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9107 -0.8800  0.0804  0.5946  4.0239 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.976e-16  5.892e-02   0.000    1.000
TDDS_p_z    1.440e-02  5.902e-02   0.244    0.807

Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared:  0.0002075, Adjusted R-squared:  -0.003276 
F-statistic: 0.05956 on 1 and 287 DF,  p-value: 0.8074

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .014, t(287) = .244, p = .807), just as with the correlation above.

  • Although this precludes the possibility that out-group avoidance mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_z, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Out-Group Avoidance (GENE_p)",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M2, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2 <- Boot(H2_M2, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M2) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE   bootMed
(Intercept) -8.8123e-16 0.00268445 0.059258 0.0010875
GENE_p_z     6.7778e-02 0.00038975 0.059761 0.0678462

confint(boot_H2_M2) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.11359766 0.1204258
GENE_p_z    -0.04811381 0.1832814

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religiosity.

• Based on the estimates, the standardized coefficient for out-group avoidance is around .068.

  • However, this relationship is not significant, as the BCI contains zero (upper = .19, lower = -.04).

Full Mediation Model

The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H2_full <- lm(CRS_z ~ TDDS_p_z + GENE_p_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_full, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_full))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H2_full))
   qqline(residuals(H2_full), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H2_full)

   TDDS_p_z GENE_p_z 
   1.000208 1.000208 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full <- Boot(H2_full, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_full) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE   bootMed
(Intercept) -8.6771e-16 0.00171383 0.057648 0.0032765
TDDS_p_z     2.4206e-01 0.00079979 0.058930 0.2414062
GENE_p_z     6.4291e-02 0.00032819 0.056867 0.0657697

confint(boot_H2_full) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.11852205 0.1083173
TDDS_p_z     0.13062951 0.3574411
GENE_p_z    -0.05479052 0.1697520

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results <- mediate(model.m = H2_M1,
        model.y = H2_full,
        treat = "TDDS_p_z",
        mediator = "GENE_p_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H2_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.000926    -0.008506         0.01    0.85    
ADE             0.242061     0.130624         0.35  <2e-16 ***
Total Effect    0.242987     0.130352         0.35  <2e-16 ***
Prop. Mediated  0.003811    -0.046003         0.06    0.85    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.

  • The coefficients for pathogen disgust (.24) and out-group avoidance (.06) are virtually identical in this model to the models above where they are stand-alone predictors. This suggests that the unique effect of pathogen disgust (and not out-group avoidance) on religiosity over-and-above each other are completely independent.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .001, BCI = [-.009, .01]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance.

Hypothesis 3: Sexual Strategies

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religiosity.

# Fitting the models
H3_M1 <- lm(SOI_a_r_z ~ TDDS_p_z, data = data)
H3_M2 <- lm(CRS_z ~ SOI_a_r_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_M1, which = 1) # For model 1

   

   plot(H3_M2, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_M1)) # For model 1

   

   hist(residuals(H3_M2)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H3_M1))
   qqline(residuals(H3_M1), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H3_M2))
   qqline(residuals(H3_M2), col = "red")

   

   • Both Models:

     • The QQ-plots for both models do not indicate that the residuals are normally distributed.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$SOI_a_r_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Restricted Sociosexual Attitudes",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M1, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1 <- Boot(H3_M1, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M1) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE    bootMed
(Intercept) -1.5075e-15 -0.0016338 0.057737 -0.0025904
TDDS_p_z     2.2087e-01 -0.0029972 0.055845  0.2187473

confint(boot_H3_M1) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1159650 0.1223375
TDDS_p_z     0.1132985 0.3316029

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)

• We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.

• Based on the estimates, the standardized coefficient for pathogen disgust is around .22.

  • In addition, this relationahip is significant, as the BCI excludes zero (upper = .33, lower = .10).

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Restricted Sociosexual Attitudes",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M2, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2 <- Boot(H3_M2, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M2) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original  bootBias   bootSE   bootMed
(Intercept) -2.4601e-16 0.0025578 0.053518 0.0036932
SOI_a_r_z    4.0934e-01 0.0016313 0.054043 0.4116882

confint(boot_H3_M2) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %     97.5 %
(Intercept) -0.1133516 0.09850889
SOI_a_r_z    0.2938363 0.50917573

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.

• Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .41.

  • In addition, this relationahip is significant, as the BCI excludes zero (upper = .51, lower = .30).

Full Mediation Model

The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.

# Fitting the full mediation model
H3_full <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_full, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_full))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H3_full))
   qqline(residuals(H3_full), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H3_full)

    TDDS_p_z SOI_a_r_z 
    1.051286  1.051286 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full <- Boot(H3_full, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_full) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original  bootBias   bootSE   bootMed
(Intercept) -2.9133e-16 0.0017123 0.053076 0.0027606
TDDS_p_z     1.6040e-01 0.0021068 0.057003 0.1635359
SOI_a_r_z    3.7391e-01 0.0014924 0.056727 0.3760744

confint(boot_H3_full) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %     97.5 %
(Intercept) -0.10951763 0.09950991
TDDS_p_z     0.05084147 0.26724821
SOI_a_r_z    0.25700020 0.47854943

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results <- mediate(model.m = H3_M1,
        model.y = H3_full,
        treat = "TDDS_p_z",
        mediator = "SOI_a_r_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H3_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             0.0826       0.0394         0.13  <2e-16 ***
ADE              0.1604       0.0515         0.26   0.008 ** 
Total Effect     0.2430       0.1304         0.35  <2e-16 ***
Prop. Mediated   0.3399       0.1682         0.67  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, both pathogen disgust and restricted sociosexual attitudes are significantly positively related to religiosity, with the unique effect of restricted sociosexual attitudes being substantially larger (B = .374, BCI =[.255, .484]) than the unique effect of pathogen disgust (B = .160, BCI = [.055, .272]).

• Mediation Analysis Results:

  • The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .08, BCI = [.038, .13]).

    • Given this result, we have evidence for the weak version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is partially mediated by a monogamous mating strategy.

  • The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .34. That is, 34% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.164, .63]).

    • Given this result, we do not have evidence for the strong version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is not fully mediated by a monogamous mating strategy.

Exploratory Analyses

Hypothesis 1: Controlling for Sex

Preliminary Simple Linear Regression Models

We will now re-test Hypothesis 1 while controlling for sex.

# Fitting the models
H1_M1_s <- lm(CONV_f_z ~ TDDS_p_z + sex, data = data)
H1_M2_s <- lm(CRS_z ~ CONV_f_z + sex, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_M1_s, which = 1) # For model 1

   

   plot(H1_M2_s, which = 1) # For model 2

   

   • The residuals vs. fitted plot for the first model is quite strange, and it suggests that there are strongly different predictions for males and females. I will differentiate the sexes by making their points different colors, to see if this is the case.
   • The residuals vs. fitted plot for the second model has a pretty clear trend towards more negative residuals at higher fitted values. We should conduct bootstrapped resampling inferences for these models.

   # Extract residuals and fitted values for model 1
   residuals_m1 <- residuals(H1_M1_s)
   fitted_m1 <- fitted(H1_M1_s)
   
   # Plot residuals vs fitted values for model 1 with different colors for sex
   plot(fitted_m1, residuals_m1, col = as.numeric(data$sex), 
        main = "Residuals vs Fitted Values (Model 1)",
        xlab = "Fitted Values", ylab = "Residuals")
   abline(h = 0, col = "red")
   
   # Add a legend to the plot
   legend("topright", legend = levels(data$sex), col = 1:length(levels(data$sex)), pch = 1)

   

   # Drop the residuals and fitted values from the environment
   rm(residuals_m1, fitted_m1)

   • Indeed, the different fitted values represent different predictions for each sex.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_M1_s)) # For model 1

   

   hist(residuals(H1_M2_s)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H1_M1_s))
   qqline(residuals(H1_M1_s), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H1_M2_s))
   qqline(residuals(H1_M2_s), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.

   • We will conduct a bootstrapped tests of the predictors for both models.

3. Multicollinearity: Variance inflation factor (VIF) for both models

   # Calculating the VIF for both models
   vif(H1_M1_s)

   TDDS_p_z      sex 
   1.038981 1.038981 

   vif(H1_M2_s)

   CONV_f_z      sex 
   1.004388 1.004388 

   • The VIF values are good here.

Summarizing the Models

Model 1

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M1_s <- Boot(H1_M1_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_M1_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
             original    bootBias   bootSE   bootMed
(Intercept)  0.068890 -0.00144952 0.087173  0.071407
TDDS_p_z     0.018766  0.00419217 0.066166  0.021675
sexFemale   -0.139225  0.00018887 0.121467 -0.139117

confint(boot_H1_M1_s) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %     97.5 %
(Intercept) -0.1101433 0.23420257
TDDS_p_z    -0.1080678 0.14762599
sexFemale   -0.3726060 0.08793543

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, TDDS_p_z, sex)

• Based on the estimate and BCI for pathogen disgust (B = .018, BCI = [-.114, .153]), controlling for sex does not change whether there is a positive relationship between pathogen disgust and traditionalism.

  • Although this precludes the possibility that adherence to tradition mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_s <- Boot(H1_M2_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_M2_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
            original   bootBias   bootSE  bootMed
(Intercept) -0.18812  0.0034793 0.075256 -0.18351
CONV_f_z     0.32702  0.0035275 0.051563  0.33122
sexFemale    0.38019 -0.0022518 0.109868  0.37625

confint(boot_H1_M2_s) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %      97.5 %
(Intercept) -0.3465750 -0.04560402
CONV_f_z     0.2115510  0.42600748
sexFemale    0.1722696  0.59940976

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, sex)

• Based on the estimates, traditionalism is significantly related to religiosity independent of sex (B = .327, BCI = [.22, .43]), and the estimate is very similar to the previous estimate without controlling for sex, suggesting that the effects are primarily independent.

Full Mediation Model

The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H1_full_s <- lm(CRS_z ~ TDDS_p_z + CONV_f_z + sex, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_full_s, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This would be problematic for parametric tests.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_full_s))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H1_full_s))
   qqline(residuals(H1_full_s), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at just around 1 SD.

   • We will conduct bootstrapping for confidence intervals to make inferences.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H1_full_s)

   TDDS_p_z CONV_f_z      sex 
   1.039334 1.004730 1.043866 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_s <- Boot(H1_full_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_full_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
            original    bootBias   bootSE  bootMed
(Intercept) -0.14722  0.00189747 0.074983 -0.14490
TDDS_p_z     0.21242 -0.00055277 0.057913  0.21240
CONV_f_z     0.32316  0.00279282 0.049039  0.32841
sexFemale    0.29753 -0.00086310 0.108512  0.29387

confint(boot_H1_full_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %       97.5 %
(Intercept) -0.29755901 -0.005976442
TDDS_p_z     0.09737364  0.322631503
CONV_f_z     0.21168286  0.409014822
sexFemale    0.09173913  0.515949990

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z, sex)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_s_mediate_results <- mediate(model.m = H1_M1_s,
        model.y = H1_full_s,
        treat = "TDDS_p_z",
        mediator = "CONV_f_z",
        covariates = "sex",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H1_s_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

(Inference Conditional on the Covariate Values Specified in `covariates')

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.00606     -0.03198         0.05   0.790    
ADE             0.21242      0.09856         0.32  <2e-16 ***
Total Effect    0.21849      0.10509         0.34   0.002 ** 
Prop. Mediated  0.02776     -0.20174         0.25   0.788    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, all three variables are significantly positively related to religiosity (with the dummy coding of sex being 1 = Female).

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = .006, BCI = [-.032, .05]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, even when controlling for sex.

Hypothesis 2: Controlling for Sex

Preliminary Simple Linear Regression Models

We will now re-test Hypothesis 2 while controlling for sex.

# Fitting the models
H2_M1_s <- lm(GENE_p_z ~ TDDS_p_z + sex, data = data)
H2_M2_s <- lm(CRS_z ~ GENE_p_z + sex, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_M1_s, which = 1) # For model 1

   

   plot(H2_M2_s, which = 1) # For model 2

   

   • The residuals vs. fitted plot for the first model seems to be partitioned by sex like the previous plot, and it suggests that there are strongly different predictions for males and females. I will differentiate the sexes by making their points different colors, to see if this is the case.
   • The second plot looks quite good.

   # Extract residuals and fitted values for model 1
   residuals_m1 <- residuals(H2_M1_s)
   fitted_m1 <- fitted(H2_M1_s)
   
   # Plot residuals vs fitted values for model 1 with different colors for sex
   plot(fitted_m1, residuals_m1, col = as.numeric(data$sex), 
        main = "Residuals vs Fitted Values (Model 1)",
        xlab = "Fitted Values", ylab = "Residuals")
   abline(h = 0, col = "red")
   
   # Add a legend to the plot
   legend("topright", legend = levels(data$sex), col = 1:length(levels(data$sex)), pch = 1)

   

   # Drop the residuals and fitted values from the environment
   rm(residuals_m1, fitted_m1)

   • Indeed, the different fitted values represent different predictions for each sex.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_M1_s)) # For model 1

   

   hist(residuals(H2_M2_s)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H2_M1_s))
   qqline(residuals(H2_M1_s), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H2_M2_s))
   qqline(residuals(H2_M2_s), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just below 2 SD.

   • We will conduct a bootstrapped tests of the predictors for both models.

3. Multicollinearity: Variance inflation factor (VIF) for both models

   # Calculating the VIF for both models
   vif(H2_M1_s)

   TDDS_p_z      sex 
   1.038981 1.038981 

   vif(H2_M2_s)

   GENE_p_z      sex 
   1.013159 1.013159 

   • The VIF values are good here.

Summarizing the Models

Model 1

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M1_s <- Boot(H2_M1_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M1_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
             original    bootBias   bootSE   bootMed
(Intercept)  0.119846  0.00029527 0.084531  0.117804
TDDS_p_z     0.037901 -0.00086675 0.060272  0.035864
sexFemale   -0.242207 -0.00040661 0.117048 -0.242947

confint(boot_H2_M1_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %      97.5 %
(Intercept) -0.03909517  0.29625715
TDDS_p_z    -0.08327969  0.15957373
sexFemale   -0.46669146 -0.01264589

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, TDDS_p_z, sex)

• Based on the estimate and BCI for pathogen disgust (B = .038, BCI = [-.084, .157]), while controlling for sex there is not a relationship between pathogen disgust and out-group avoidance.

  • Although this precludes the possibility that out-group avoidance mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_s <- Boot(H2_M2_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M2_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
            original    bootBias   bootSE   bootMed
(Intercept) -0.17669  0.00383719 0.080908 -0.174657
GENE_p_z     0.08816  0.00058913 0.058323  0.089972
sexFemale    0.35710 -0.00291569 0.118240  0.350989

confint(boot_H2_M2_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %      97.5 %
(Intercept) -0.33163460 -0.01089098
GENE_p_z    -0.03812488  0.19549662
sexFemale    0.13230476  0.58667840

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, sex)

• Based on the estimates, out-group avoidance is not significantly related to religiosity independent of sex (B = .088, BCI = [-.025, .214]).

Full Mediation Model

The fact that pathogen disgust does not predict out-group avoidance and that out-group avoidance does not predict religiosity indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H2_full_s <- lm(CRS_z ~ TDDS_p_z + GENE_p_z + sex, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_full_s, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This would be problematic for parametric tests.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_full_s))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H2_full_s))
   qqline(residuals(H2_full_s), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at just above -1 SD and just above 1 SD.

   • We will conduct bootstrapping for confidence intervals to make inferences.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H2_full_s)

   TDDS_p_z GENE_p_z      sex 
   1.040438 1.014580 1.053911 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_s <- Boot(H2_full_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_full_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
             original    bootBias   bootSE   bootMed
(Intercept) -0.134569  0.00215572 0.080493 -0.131898
TDDS_p_z     0.215448  0.00086485 0.061324  0.214958
GENE_p_z     0.080198  0.00054884 0.056106  0.080718
sexFemale    0.271962 -0.00137026 0.117763  0.263523

confint(boot_H2_full_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.29354508 0.0244336
TDDS_p_z     0.09706498 0.3370680
GENE_p_z    -0.03259225 0.1865872
sexFemale    0.05112016 0.5181156

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z, sex)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_s_mediate_results <- mediate(model.m = H2_M1_s,
        model.y = H2_full_s,
        treat = "TDDS_p_z",
        mediator = "GENE_p_z",
        covariates = "sex",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H2_s_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

(Inference Conditional on the Covariate Values Specified in `covariates')

               Estimate 95% CI Lower 95% CI Upper p-value   
ACME            0.00304     -0.00774         0.02   0.588   
ADE             0.21545      0.10085         0.33   0.002 **
Total Effect    0.21849      0.10509         0.34   0.002 **
Prop. Mediated  0.01391     -0.04180         0.10   0.586   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, both pathogen disgust and sex, but not out-group avoidance, are significantly positively related to religiosity (with the dummy coding of sex being 1 = Female).

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .003, BCI = [-.007, .02]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even when controlling for sex.

Hypothesis 3: Controlling for Sex

Preliminary Simple Linear Regression Models

We will now re-test Hypothesis 3 while controlling for sex.

# Fitting the models
H3_M1_s <- lm(SOI_a_r_z ~ TDDS_p_z + sex, data = data)
H3_M2_s <- lm(CRS_z ~ SOI_a_r_z + sex, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_M1_s, which = 1) # For model 1

   

   plot(H3_M2_s, which = 1) # For model 2

   

   • Unlike the previous two plots, the predicted values for model 1 are not very different by sex (or at least not obviously). There is a bit of downward trend, so we will bootstrap for inferences of the coefficients.
   • The second plot is quite similar to the first.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_M1_s)) # For model 1

   

   hist(residuals(H3_M2_s)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H3_M1_s))
   qqline(residuals(H3_M1_s), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H3_M2_s))
   qqline(residuals(H3_M2_s), col = "red")

   

   • Both Models:

     • The QQ-plots for both models are not as bad as in previous models, but there is some deviation from normality towards the ends of the distributions.

   • We will conduct a bootstrapped tests of the predictors for both models.

3. Multicollinearity: Variance inflation factor (VIF) for both models

   # Calculating the VIF for both models
   vif(H3_M1_s)

   TDDS_p_z      sex 
   1.038981 1.038981 

   vif(H3_M2_s)

   SOI_a_r_z       sex 
    1.038623  1.038623 

   • The VIF values are good here.

Summarizing the Models

Model 1

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_s <- Boot(H3_M1_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M1_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
            original    bootBias   bootSE  bootMed
(Intercept) -0.15403 -3.7669e-05 0.080601 -0.15429
TDDS_p_z     0.19067 -2.4728e-03 0.057487  0.18827
sexFemale    0.31129 -3.3652e-03 0.115848  0.30770

confint(boot_H3_M1_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %      97.5 %
(Intercept) -0.30713596 0.004964965
TDDS_p_z     0.07973518 0.305882636
sexFemale    0.09905334 0.546865791

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z, sex)

• Based on the estimate and BCI for pathogen disgust (B = .191, BCI = [.076, .309]) and sex (B = .311, BCI = [.079, .541]), there is a significant relationship between both variables and restricted sociosexual attitudes.

Model 2

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_s <- Boot(H3_M2_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M2_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
             original   bootBias   bootSE  bootMed
(Intercept) -0.092212  0.0034832 0.078908 -0.08793
SOI_a_r_z    0.391337  0.0023672 0.054617  0.39386
sexFemale    0.186359 -0.0020876 0.114311  0.18463

confint(boot_H3_M2_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %     97.5 %
(Intercept) -0.24883828 0.06126489
SOI_a_r_z    0.27960214 0.49433536
sexFemale   -0.02594773 0.40381904

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, sex)

• Based on the estimates, restricted sociosexual attitudes is significantly related to religiosity independent of sex (B = .391, BCI = [.271, .493), and sex drops just below significance as a predictor of religiosity (B = .186, BCI = [-.026, .417]).

Full Mediation Model

# Fitting the full mediation model
H3_full_s <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z + sex, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_full_s, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This would be problematic for parametric tests.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_full_s))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H3_full_s))
   qqline(residuals(H3_full_s), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not quite normally distributed.

   • We will conduct bootstrapping for confidence intervals to make inferences.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H3_full_s)

    TDDS_p_z SOI_a_r_z       sex 
    1.078165  1.077794  1.065178 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_s <- Boot(H3_full_s, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_full_s) # Getting the coefficients

Number of bootstrap replications R = 1000 
             original   bootBias   bootSE   bootMed
(Intercept) -0.069064  0.0022568 0.078290 -0.066787
TDDS_p_z     0.149296  0.0019682 0.058524  0.152164
SOI_a_r_z    0.362880  0.0020451 0.057094  0.365231
sexFemale    0.139577 -0.0013943 0.113992  0.137545

confint(boot_H3_full_s) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %     97.5 %
(Intercept) -0.22670927 0.07918258
TDDS_p_z     0.03465778 0.25852694
SOI_a_r_z    0.25101165 0.46914310
sexFemale   -0.07659668 0.37316434

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z, sex)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_s_mediate_results <- mediate(model.m = H3_M1_s,
        model.y = H3_full_s,
        treat = "TDDS_p_z",
        mediator = "SOI_a_r_z",
        covariates = "sex",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H3_s_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

(Inference Conditional on the Covariate Values Specified in `covariates')

               Estimate 95% CI Lower 95% CI Upper p-value   
ACME             0.0692       0.0290         0.12   0.002 **
ADE              0.1493       0.0366         0.26   0.012 * 
Total Effect     0.2185       0.1051         0.34   0.002 **
Prop. Mediated   0.3167       0.1444         0.69   0.004 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, both pathogen disgust and restricted sociosexual attitudes, but not sex, are significantly positively related to religiosity (with the dummy coding of sex being 1 = Female).

• Mediation Analysis Results:

• The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .069, BCI = [.027, .12]).

  • Given this result, we, again, have evidence for the weak version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is partially mediated by a monogamous mating strategy, even while controlling for sex.

• The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .32. That is, 32% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.136, .62]).

  • Given this result, we do not have evidence for the strong version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is not fully mediated by a monogamous mating strategy.

Hypothesis 1: Religious Importance

Preliminary Simple Linear Regression Models

Now we will retest each hypothesis with religious importance as the outcome variable, as a check for robustness.

As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religious importance.

# Fitting the models
H1_M1_RI <- lm(CONV_f_z ~ TDDS_p_z, data = data)
H1_M2_RI <- lm(religious_importance_z ~ CONV_f_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_M1_RI, which = 1) # For model 1

   

   plot(H1_M2_RI, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the first model, but there is a trend towards more negative residual values at the higher fitted values of the model for model 2. This is problematic for parametric inference.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_M1_RI)) # For model 1

   

   hist(residuals(H1_M2_RI)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H1_M1_RI))
   qqline(residuals(H1_M1_RI), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H1_M2_RI))
   qqline(residuals(H1_M2_RI), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CONV_f_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Traditionalism",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M1_RI, col = "black", lwd = 1)

# Summarizing the model
summary(H1_M1_RI)

Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4517 -0.6429  0.0402  0.6287  2.8524 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.379e-17  5.893e-02   0.000    1.000
TDDS_p_z     5.260e-03  5.903e-02   0.089    0.929

Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared:  2.766e-05, Adjusted R-squared:  -0.003457 
F-statistic: 0.00794 on 1 and 287 DF,  p-value: 0.9291

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = .005, t(287) = .089, p = .929), just as with the correlation above.

  • Although this precludes the possibility that adherence to tradition mediates the relationship between disgust and religious importance, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$CONV_f_z, data$religious_importance_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Traditionalism",
     ylab = "Religious Importance",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M2_RI, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
religious_importance_z <- as.numeric(data$religious_importance_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_RI <- Boot(H1_M2_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_M2_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original  bootBias   bootSE   bootMed
(Intercept) 1.7368e-16 0.0029611 0.055114 0.0011843
CONV_f_z    3.3199e-01 0.0044058 0.052992 0.3353958

confint(boot_H1_M2_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1068043 0.1145093
CONV_f_z     0.2189193 0.4346032

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, religious_importance_z)

• We can see from the Scatter Plot that there is a positive relationship between traditionalism and religious importance.

• Based on the estimates, the standardized coefficient for traditionalism is around .33.

  • In addition, this relationship is significant, as the BCIs do not contain zero (upper = .44, lower = .22).

Full Mediation Model

The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religious importance, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H1_full_RI <- lm(religious_importance_z ~ TDDS_p_z + CONV_f_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_full_RI, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This is problematic for parametric inference.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_full_RI))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H1_full_RI))
   qqline(residuals(H1_full_RI), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H1_full_RI)

   TDDS_p_z CONV_f_z 
   1.000028 1.000028 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_RI <- Boot(H1_full_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_full_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE  bootMed
(Intercept) 1.8511e-16  0.00205054 0.053908 0.002309
TDDS_p_z    2.1682e-01 -0.00097327 0.057490 0.214199
CONV_f_z    3.3085e-01  0.00389128 0.050288 0.334003

confint(boot_H1_full_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1082920 0.1108538
TDDS_p_z     0.1105920 0.3414131
CONV_f_z     0.2227669 0.4216086

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, religious_importance_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_RI_mediate_results <- mediate(model.m = H1_M1_RI,
        model.y = H1_full_RI,
        treat = "TDDS_p_z",
        mediator = "CONV_f_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H1_RI_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.00174     -0.03650         0.05    0.96    
ADE             0.21682      0.09311         0.33  <2e-16 ***
Total Effect    0.21856      0.10172         0.34  <2e-16 ***
Prop. Mediated  0.00796     -0.22524         0.24    0.96    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religious importance.

  • The coefficients for pathogen disgust (.22) and traditionalism (.33) are similar to the model with the CRS.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religious importance through traditionalism (B = .002, BCI = [-.038, .04]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religious importance is mediated by traditionalism.

Hypothesis 2: Religious Importance

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religious importance.

# Fitting the models
H2_M1_RI <- lm(GENE_p_z ~ TDDS_p_z, data = data)
H2_M2_RI <- lm(religious_importance_z ~ GENE_p_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_M1_RI, which = 1) # For model 1

   

   plot(H2_M2_RI, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_M1_RI)) # For model 1

   

   hist(residuals(H2_M2_RI)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H2_M1_RI))
   qqline(residuals(H2_M1_RI), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H2_M2_RI))
   qqline(residuals(H2_M2_RI), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are relatively normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Out-Group Avoidance (GENE_p)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M1_RI, col = "black", lwd = 1)

# Summarizing the model
summary(H2_M1_RI)

Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9107 -0.8800  0.0804  0.5946  4.0239 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.976e-16  5.892e-02   0.000    1.000
TDDS_p_z    1.440e-02  5.902e-02   0.244    0.807

Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared:  0.0002075, Adjusted R-squared:  -0.003276 
F-statistic: 0.05956 on 1 and 287 DF,  p-value: 0.8074

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .014, t(287) = .244, p = .807)—same as the model above.

  • Although this precludes the possibility that out-group avoidance mediates the relationship between disgust and religious importance, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_z, data$religious_importance_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Out-Group Avoidance (GENE_p)",
     ylab = "Religious Importance",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M2_RI, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
religious_importance_z <- as.numeric(data$religious_importance_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_RI <- Boot(H2_M2_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M2_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE   bootMed
(Intercept) 1.4098e-16  0.00271370 0.059242 0.0014941
GENE_p_z    7.4934e-02 -0.00037603 0.060076 0.0762719

confint(boot_H2_M2_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.10931327 0.1249488
GENE_p_z    -0.04800801 0.1879390

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, religious_importance_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religious importance.

• Based on the estimates, the standardized coefficient for traditionalism is around .07.

  • However, this relationship is not significant, as the BCI contains zero (upper = .20, lower = -.04).

Full Mediation Model

The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religious importance indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religious importance, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H2_full_RI <- lm(religious_importance_z ~ TDDS_p_z + GENE_p_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_full_RI, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_full_RI))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H2_full_RI))
   qqline(residuals(H2_full_RI), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H2_full_RI)

   TDDS_p_z GENE_p_z 
   1.000208 1.000208 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_RI <- Boot(H2_full_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_full_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE   bootMed
(Intercept) 1.5313e-16  0.00183362 0.057942 0.0011976
TDDS_p_z    2.1752e-01  0.00041664 0.059876 0.2144363
GENE_p_z    7.1801e-02 -0.00026517 0.057177 0.0725622

confint(boot_H2_full_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.10682993 0.1247495
TDDS_p_z     0.10889408 0.3438263
GENE_p_z    -0.04473746 0.1872367

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, religious_importance_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_RI_mediate_results <- mediate(model.m = H2_M1_RI,
        model.y = H2_full_RI,
        treat = "TDDS_p_z",
        mediator = "GENE_p_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H2_RI_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.00103     -0.00929         0.02    0.84    
ADE             0.21752      0.10126         0.34  <2e-16 ***
Total Effect    0.21856      0.10172         0.34  <2e-16 ***
Prop. Mediated  0.00473     -0.04825         0.08    0.84    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religious importance.

  • The coefficients for pathogen disgust (.22) and out-group avoidance (.07) are virtually identical in the model with CRS.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religious importance through out-group avoidance (B = .001, BCI = [-.009, .01]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religious importance is mediated by out-group avoidance.

Hypothesis 3: Religious Importance

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religious importance.

# Fitting the models
H3_M1_RI <- lm(SOI_a_r_z ~ TDDS_p_z, data = data)
H3_M2_RI <- lm(religious_importance_z ~ SOI_a_r_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_M1_RI, which = 1) # For model 1

   

   plot(H3_M2_RI, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for model 1, but model two has a pretty decent trend toward more negative residual values at higher fitted values.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_M1_RI)) # For model 1

   

   hist(residuals(H3_M2_RI)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H3_M1_RI))
   qqline(residuals(H3_M1_RI), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H3_M2_RI))
   qqline(residuals(H3_M2_RI), col = "red")

   

   • Both Models:

     • The QQ-plots for both models do not indicate that the residuals are normally distributed.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$SOI_a_r_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Restricted Sociosexual Attitudes",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M1_RI, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_RI <- Boot(H3_M1_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M1_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE    bootMed
(Intercept) -1.5075e-15 -0.0016338 0.057737 -0.0025904
TDDS_p_z     2.2087e-01 -0.0029972 0.055845  0.2187473

confint(boot_H3_M1_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1159650 0.1223375
TDDS_p_z     0.1132985 0.3316029

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)

• We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.

• Based on the estimates, the standardized coefficient for pathogen disgust is around .22.

  • In addition, this relationship is significant, as the BCI excludes zero (upper = .33, lower = .10).

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$religious_importance_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Restricted Sociosexual Attitudes",
     ylab = "Religious Importance",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M2_RI, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
religious_importance_z <- as.numeric(data$religious_importance_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_RI <- Boot(H3_M2_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M2_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original  bootBias   bootSE   bootMed
(Intercept) 7.4203e-16 0.0024378 0.054107 0.0011108
SOI_a_r_z   3.8592e-01 0.0025175 0.055091 0.3877059

confint(boot_H3_M2_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1054531 0.1082780
SOI_a_r_z    0.2755024 0.4862312

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, religious_importance_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religious importance.

• Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .39.

  • In addition, this relationship is significant, as the BCI excludes zero (upper = .49, lower = .27).

Full Mediation Model

The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religious importance indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religious importance. Now we will run the mediation analysis.

# Fitting the full mediation model
H3_full_RI <- lm(religious_importance_z ~ TDDS_p_z + SOI_a_r_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_full_RI, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_full_RI))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H3_full_RI))
   qqline(residuals(H3_full_RI), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H3_full_RI)

    TDDS_p_z SOI_a_r_z 
    1.051286  1.051286 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_RI <- Boot(H3_full_RI, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_full_RI) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original  bootBias   bootSE   bootMed
(Intercept) 7.0243e-16 0.0017189 0.053773 0.0011696
TDDS_p_z    1.4015e-01 0.0018125 0.057063 0.1415711
SOI_a_r_z   3.5496e-01 0.0021128 0.057777 0.3581158

confint(boot_H3_full_RI) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1063690 0.1042462
TDDS_p_z     0.0341517 0.2534335
SOI_a_r_z    0.2322735 0.4558705

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, religious_importance_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_RI_mediate_results <- mediate(model.m = H3_M1_RI,
        model.y = H3_full_RI,
        treat = "TDDS_p_z",
        mediator = "SOI_a_r_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H3_RI_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             0.0784       0.0371         0.13  <2e-16 ***
ADE              0.1402       0.0255         0.25   0.016 *  
Total Effect     0.2186       0.1017         0.34  <2e-16 ***
Prop. Mediated   0.3587       0.1815         0.78  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, both pathogen disgust and restricted sociosexual attitudes are significantly positively related to religious importance, with the unique effect of restricted sociosexual attitudes being substantially larger (B = .355, BCI =[.237, .467]) than the unique effect of pathogen disgust (B = .140, BCI = [.036, .253]).

• Mediation Analysis Results:

  • The ACME indicates that there is a significant indirect effect of pathogen disgust on religious importance through a monogamous mating strategy (B = .078, BCI = [.037, .13]).

    • Given this result, we have evidence for the weak version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is partially mediated by a monogamous mating strategy.

  • The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .36. That is, 36% of the covariance between pathogen disgust and religious importance seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.176, .71]).

    • Given this result, we do not have evidence for the strong version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is not fully mediated by a monogamous mating strategy.

Hypothesis 2: Removing GENE Preferences Item 1

In the Internal Reliability analysis of the Preferences Subscale of SFGENE-7, the drop alpha statistics indicated that removing Item 1 from the subscale would improve the internal reliability from an alpha of around .5 to around .8. Given that a lack of internal reliability may attenuate relationships between variables in correlational analyses, I will now conduct the same analysis for Hypothesis 2 but with Item 1 removed.

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_2_z) and that out-group avoidance influences religiosity.

# Fitting the models
H2_M1_2 <- lm(GENE_p_2_z ~ TDDS_p_z, data = data)
H2_M2_2 <- lm(CRS_z ~ GENE_p_2_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_M1_2, which = 1) # For model 1

   

   plot(H2_M2_2, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_M1_2)) # For model 1

   

   hist(residuals(H2_M2_2)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H2_M1_2))
   qqline(residuals(H2_M1_2), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H2_M2_2))
   qqline(residuals(H2_M2_2), col = "red")

   

   • Both Models:

     • The histograms and QQ-plots indicate that for both models the residuals are not normally distributed, although in different ways.

     • We will need to bootstrapping for testing the predictors for these models.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_2_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Out-Group Avoidance (GENE_p without Item 1)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M1_2, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M1_2 <- Boot(H2_M1_2, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M1_2) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE    bootMed
(Intercept) -1.2029e-16 -0.00018739 0.059877 -0.0018552
TDDS_p_z    -2.5360e-02  0.00068765 0.059533 -0.0261568

confint(boot_H2_M1_2) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %     97.5 %
(Intercept) -0.1102701 0.12771218
TDDS_p_z    -0.1342644 0.09246806

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, TDDS_p_z)

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = -.025, BCI = [-.133, .106])—just as with the analysis above.

  • Although this precludes the possibility that out-group avoidance mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_2_z, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Out-Group Avoidance (GENE_p without Item 1)",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M2_2, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_2 <- Boot(H2_M2_2, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M2_2) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE    bootMed
(Intercept) -8.6789e-16  2.2888e-03 0.059385 0.00088737
GENE_p_2_z   1.3388e-05 -2.0226e-05 0.062432 0.00152762

confint(boot_H2_M2_2) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1153198 0.1206342
GENE_p_2_z  -0.1212795 0.1206176

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, CRS_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religiosity.

• Based on the estimates, the standardized coefficient for out-group avoidance is around B = .00001 (BCI = [-.119, .123]).

Full Mediation Model

The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H2_2_full <- lm(CRS_z ~ TDDS_p_z + GENE_p_2_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_2_full, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_2_full))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H2_2_full))
   qqline(residuals(H2_2_full), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H2_2_full)

     TDDS_p_z GENE_p_2_z 
     1.000644   1.000644 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_2_full <- Boot(H2_2_full, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_2_full) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE   bootMed
(Intercept) -8.5427e-16  0.00137003 0.057776 0.0028594
TDDS_p_z     2.4314e-01  0.00028211 0.059630 0.2416846
GENE_p_2_z   6.1794e-03 -0.00026584 0.057790 0.0055259

confint(boot_H2_2_full) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1188244 0.1086140
TDDS_p_z     0.1337811 0.3610994
GENE_p_2_z  -0.1023601 0.1179598

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_2_mediate_results <- mediate(model.m = H2_M1_2,
        model.y = H2_2_full,
        treat = "TDDS_p_z",
        mediator = "GENE_p_2_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H2_2_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.000157    -0.006777         0.01    0.92    
ADE             0.243143     0.129699         0.35  <2e-16 ***
Total Effect    0.242987     0.130352         0.35  <2e-16 ***
Prop. Mediated -0.000645    -0.034402         0.05    0.92    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.

  • The coefficients for pathogen disgust (.24) and out-group avoidance (.006) are similar to their values with Item 1 included.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = -.0002, BCI = [-.006, .01]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance.

Data Preparation for Standardized Residuals

Now we will rerun the analyses after ensuring that covariance between our mediating variables could not account for the results above. That is, we will regress each of the mediator on the other two and save standardized residuals. Then, we will run the mediation analyses again. This will help to explore whether any mediation effects (or null effects) in our previous analyses could have been due to covariation between our mediating variables. First, let us calculate our standardized residuals.

# Create the regression models for standardized residuals
  # For traditionalism
res_model_CONV_f <- lm(CONV_f ~ SOI_a_r + GENE_p, data = data)
  # For out-group avoidance
res_model_GENE_p <- lm(GENE_p ~ CONV_f + SOI_a_r, data = data)
  # For sexual strategies
res_model_SOI_a_r <- lm(SOI_a_r ~ GENE_p + CONV_f, data = data)

# Add the standardized residuals to the new data frame, aligning them with the full rows
  # Initializing new variables with NA values first
data$res_CONV_f <- NA   # Initialize new variable with NAs for traditionalism
data$res_GENE_p <- NA   # Initialize new variable with NAs for out-group avoidance
data$res_SOI_a_r <- NA
  # Add the standardized residuals to the res_CONV_f variable
data$res_CONV_f <- rstandard(res_model_CONV_f)
  # Add the standardized residuals to the res_GENE_p variable
data$res_GENE_p <- rstandard(res_model_GENE_p)
  # Add the standardized residuals to the res_SOI_a_r variable
data$res_SOI_a_r <- rstandard(res_model_SOI_a_r)

# Ensure that the variables are plain numeric vectors
data$res_CONV_f <- as.numeric(data$res_CONV_f)
data$res_GENE_p <- as.numeric(data$res_GENE_p)
data$res_SOI_a_r <- as.numeric(data$res_SOI_a_r)

Hypothesis 1: Standardized Residuals

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (in this case the standardized residuals version; res_CONV_f) and that adherence to traditional practices influences religiosity.

# Fitting the models
H1_M1_res <- lm(res_CONV_f ~ TDDS_p_z, data = data)
H1_M2_res <- lm(CRS_z ~ res_CONV_f, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_M1_res, which = 1) # For model 1

   

   plot(H1_M2_res, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the first model, but there is a trend towards more negative residual values at the higher fitted values of the model for model 2. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_M1_res)) # For model 1

   

   hist(residuals(H1_M2_res)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H1_M1_res))
   qqline(residuals(H1_M1_res), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H1_M2_res))
   qqline(residuals(H1_M2_res), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$res_CONV_f,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Traditionalism Standardized Residuals",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M1_res, col = "black", lwd = 1)

# Summarizing the model
summary(H1_M1_res)

Call:
lm(formula = res_CONV_f ~ TDDS_p_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5869 -0.6302  0.0261  0.6106  2.9054 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0002967  0.0590112  -0.005    0.996
TDDS_p_z    -0.0272508  0.0591136  -0.461    0.645

Residual standard error: 1.003 on 287 degrees of freedom
Multiple R-squared:  0.0007399, Adjusted R-squared:  -0.002742 
F-statistic: 0.2125 on 1 and 287 DF,  p-value: 0.6452

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = -.027, t(287) = .461, p = .645), just as with the analysis above.

  • Although this precludes the possibility that adherence to tradition mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$res_CONV_f, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Traditionalism Standardized Residuals",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M2_res, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
res_CONV_f <- as.numeric(data$res_CONV_f)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_res <- Boot(H1_M2_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_M2_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original  bootBias   bootSE   bootMed
(Intercept) 7.5754e-05 0.0026181 0.056525 0.0028825
res_CONV_f  2.5535e-01 0.0034721 0.054933 0.2591966

confint(boot_H1_M2_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1156756 0.1095724
res_CONV_f   0.1345532 0.3565031

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_CONV_f, CRS_z)

• We can see from the Scatter Plot that there is a positive relationship between traditionalism and religiosity.

• Based on the estimates, the standardized coefficient for traditionalism is around .25.

  • In addition, this relationship is significant, as the BCIs do not contain zero (upper = .36, lower = .14).

Full Mediation Model

The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H1_full_res <- lm(CRS_z ~ TDDS_p_z + res_CONV_f, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_full_res, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_full_res))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H1_full_res))
   qqline(residuals(H1_full_res), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and around 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H1_full_res)

     TDDS_p_z res_CONV_f 
      1.00074    1.00074 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
res_CONV_f <- as.numeric(data$res_CONV_f)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_res <- Boot(H1_full_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_full_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE   bootMed
(Intercept) 7.7769e-05  0.00163568 0.054903 0.0027247
TDDS_p_z    2.5013e-01 -0.00066573 0.057291 0.2487006
res_CONV_f  2.6214e-01  0.00245588 0.050525 0.2645166

confint(boot_H1_full_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1114009 0.1040761
TDDS_p_z     0.1420620 0.3667127
res_CONV_f   0.1571168 0.3577343

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_CONV_f, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results_res <- mediate(model.m = H1_M1_res,
        model.y = H1_full_res,
        treat = "TDDS_p_z",
        mediator = "res_CONV_f",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H1_mediate_results_res)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.00714     -0.03622         0.03    0.63    
ADE             0.25013      0.13001         0.36  <2e-16 ***
Total Effect    0.24299      0.13035         0.35  <2e-16 ***
Prop. Mediated -0.02940     -0.19132         0.12    0.63    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religiosity.

  • The coefficients for pathogen disgust (.25) and traditionalism (.26) are close to how they are in the models above where they are stand-alone predictors. This suggests that the unique effects of pathogen disgust and traditionalism on religiosity over-and-above each other are independent effects.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = -.007, BCI = [-.038, .03]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, even after removing covariation with the other potentially mediating variables.

Hypothesis 2: Standardized Residuals

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., the standardized residual version of res_GENE_p) and that out-group avoidance influences religiosity.

# Fitting the models
H2_M1_res <- lm(res_GENE_p ~ TDDS_p_z, data = data)
H2_M2_res <- lm(CRS_z ~ res_GENE_p, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_M1_res, which = 1) # For model 1

   

   plot(H2_M2_res, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_M1_res)) # For model 1

   

   hist(residuals(H2_M2_res)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H2_M1_res))
   qqline(residuals(H2_M1_res), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H2_M2_res))
   qqline(residuals(H2_M2_res), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are relatively normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$res_GENE_p,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Out-Group Avoidance (res_GENE_p) Standardized Residuals",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M1_res, col = "black", lwd = 1)

# Summarizing the model
summary(H2_M1_res)

Call:
lm(formula = res_GENE_p ~ TDDS_p_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1221 -0.7581 -0.0587  0.7029  3.9597 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0001159  0.0590329  -0.002    0.998
TDDS_p_z    -0.0175214  0.0591353  -0.296    0.767

Residual standard error: 1.004 on 287 degrees of freedom
Multiple R-squared:  0.0003058, Adjusted R-squared:  -0.003177 
F-statistic: 0.08779 on 1 and 287 DF,  p-value: 0.7672

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = -.017, t(287) = -.296, p = .767), just as with the correlation above.

  • Although this precludes the possibility that out-group avoidance mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$res_GENE_p, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Out-Group Avoidance (GENE_p)",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M2_res, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
res_GENE_p <- as.numeric(data$res_GENE_p)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_res <- Boot(H2_M2_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M2_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE    bootMed
(Intercept) -1.1983e-06  2.4874e-03 0.059214  0.0011912
res_GENE_p  -1.0335e-02 -9.9515e-06 0.061660 -0.0118285

confint(boot_H2_M2_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1148799 0.1187364
res_GENE_p  -0.1243771 0.1173597

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_GENE_p, CRS_z)

• We can see from the Scatter Plot that there no relationship.

• Based on the estimates, the standardized coefficient for out-group avoidance is around .01.

  • This relationship is not significant, as the BCI contains zero (upper = .13, lower = -.12).

Full Mediation Model

The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H2_full_res <- lm(CRS_z ~ TDDS_p_z + res_GENE_p, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_full_res, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_full_res))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H2_full_res))
   qqline(residuals(H2_full_res), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H2_full_res)

     TDDS_p_z res_GENE_p 
     1.000306   1.000306 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
res_GENE_p <- as.numeric(data$res_GENE_p)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_res <- Boot(H2_full_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_full_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE    bootMed
(Intercept) -7.0685e-07  0.00151630 0.057644  0.0022743
TDDS_p_z     2.4288e-01  0.00069823 0.059676  0.2411531
res_GENE_p  -6.0965e-03 -0.00012461 0.058710 -0.0056310

confint(boot_H2_full_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1181186 0.1086985
TDDS_p_z     0.1325909 0.3610056
res_GENE_p  -0.1214375 0.1088764

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_GENE_p, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results_res <- mediate(model.m = H2_M1_res,
        model.y = H2_full_res,
        treat = "TDDS_p_z",
        mediator = "res_GENE_p",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H2_mediate_results_res)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.000107    -0.006742         0.01    0.87    
ADE             0.242880     0.130981         0.36  <2e-16 ***
Total Effect    0.242987     0.130352         0.35  <2e-16 ***
Prop. Mediated  0.000440    -0.027246         0.04    0.87    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.

  • The coefficients for pathogen disgust (.24) and out-group avoidance (.-.006) are very similar in this model to the models above where they are stand-alone predictors.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .0001, BCI = [-.006, .01]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even when using the standardized residuals.

Hypothesis 3: Standardized Residuals

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; in this case the standardized residuals version res_SOI_a_r) and that a monogamous mating strategy influences religiosity.

# Fitting the models
H3_M1_res <- lm(res_SOI_a_r ~ TDDS_p_z, data = data)
H3_M2_res <- lm(CRS_z ~ res_SOI_a_r, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_M1_res, which = 1) # For model 1

   

   plot(H3_M2_res, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the first model, but there is a negative trend for residuals for the second model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_M1_res)) # For model 1

   

   hist(residuals(H3_M2_res)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H3_M1_res))
   qqline(residuals(H3_M1_res), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H3_M2_res))
   qqline(residuals(H3_M2_res), col = "red")

   

   • Both Models:

     • The QQ-plots for both models do not indicate that the residuals are normally distributed.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$res_SOI_a_r,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Restricted Sociosexual Attitudes",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M1_res, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
res_SOI_a_r <- as.numeric(data$res_SOI_a_r)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_res <- Boot(H3_M1_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M1_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original   bootBias   bootSE    bootMed
(Intercept) 4.5205e-05 -0.0015333 0.057547 -0.0027783
TDDS_p_z    2.2359e-01 -0.0035998 0.054417  0.2211374

confint(boot_H3_M1_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1134134 0.1188434
TDDS_p_z     0.1248918 0.3356576

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_SOI_a_r, TDDS_p_z)

• We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.

• Based on the estimates, the standardized coefficient for pathogen disgust is around .22.

  • In addition, this relationahip is significant, as the BCI excludes zero (upper = .34, lower = .11).

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$res_SOI_a_r, data$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Restricted Sociosexual Attitudes",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M2_res, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
res_SOI_a_r <- as.numeric(data$res_SOI_a_r)
CRS_z <- as.numeric(data$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_res <- Boot(H3_M2_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M2_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original  bootBias   bootSE   bootMed
(Intercept) -1.6429e-05 0.0024569 0.054711 0.0038558
res_SOI_a_r  3.6343e-01 0.0020824 0.056629 0.3659118

confint(boot_H3_M2_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1119157 0.1047078
res_SOI_a_r  0.2401447 0.4617742

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_SOI_a_r, CRS_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.

• Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .36.

  • In addition, this relationahip is significant, as the BCI excludes zero (upper = .47, lower = .24).

Full Mediation Model

The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.

# Fitting the full mediation model
H3_full_res <- lm(CRS_z ~ TDDS_p_z + res_SOI_a_r, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_full_res, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_full_res))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H3_full_res))
   qqline(residuals(H3_full_res), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H3_full_res)

      TDDS_p_z res_SOI_a_r 
      1.052405    1.052405 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
res_SOI_a_r <- as.numeric(data$res_SOI_a_r)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_res <- Boot(H3_full_res, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_full_res) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original  bootBias   bootSE   bootMed
(Intercept) -1.4715e-05 0.0015772 0.054194 0.0017929
TDDS_p_z     1.7020e-01 0.0020435 0.058640 0.1728617
res_SOI_a_r  3.2552e-01 0.0020325 0.059172 0.3280743

confint(boot_H3_full_res) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1088316 0.1041097
TDDS_p_z     0.0562329 0.2816977
res_SOI_a_r  0.2050904 0.4359742

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_SOI_a_r, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results_res <- mediate(model.m = H3_M1_res,
        model.y = H3_full_res,
        treat = "TDDS_p_z",
        mediator = "res_SOI_a_r",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H3_mediate_results_res)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             0.0728       0.0346         0.12  <2e-16 ***
ADE              0.1702       0.0603         0.27   0.006 ** 
Total Effect     0.2430       0.1304         0.35  <2e-16 ***
Prop. Mediated   0.2995       0.1471         0.60  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, both pathogen disgust and restricted sociosexual attitudes are significantly positively related to religiosity, with the unique effect of restricted sociosexual attitudes being substantially larger (B = .325, BCI =[.202, .441]) than the unique effect of pathogen disgust (B = .17, BCI = [.062, .289]).

• Mediation Analysis Results:

  • The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .07, BCI = [.03, .12]).

    • Given this result, we have evidence for the weak version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is partially mediated by a monogamous mating strategy.

  • The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .30. That is, 30% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.13, .57]).

    • Given this result, we do not have evidence for the strong version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is not fully mediated by a monogamous mating strategy.

Ultimately, the results of reconducting the analysis with standardized residuals reinforces the results of the analysis thus far.

Exploratory Analyses for RBB Reviewers

After receiving helpful suggestions from the reviewers at Religion, Brain and Behavior, we have decided to conduct a number of exploratory analyses that will disambiguate the findings here.

1. Reviewer 1 pointed out that almost half the sample indicate no religious affiliation and that we should present results with and without this group. This will be helpful to better characterize our results. Do individual differences in disgust sensitivity relate, through sexual strategies, to whether someone is or is not religious? Or do individual differences in disgust sensitivity relate, through sexual strategies, to someone’s degree of religiosity? By performing the analyses again on only the explicitly religious participants, we can address these questions.
2. Reviewer 1 also points out that BIS reactivity may encourage traditional cultural values (like religion) as a means of mitigating social behaviors (like sex) that are more likely to result in the spread of infectious disease. To assess this, they suggest analyzing whether sociosexual attitudes mediate the relation between sexual disgust and traditionalism, as traditionalism shows a bivariate relationship with sexual disgust above. We will address this with additional analyses.
3. Lastly, Reviewer 2 expressed concerns about population differences in disgust sensitivity and religiosity spuriously creating a relationship between disgust and religiosity at the individual level. That is, even if there is a true null relationshiop between disgust and religiosity at the individual-level within countries, the presence of one or more sampled nations where there is high disgust and high religiosity levels could result in observing a relationship between disgust and religiosity at the individual level that simply doesn’t exist. Reviewer 2 does not see a resolution to this issue. We will explore different models (e.g., a multilevel model with random intercepts to account for regional variation in the data) to try to determine whether this is the case.

Data Preparation for Religious Subgroup Analyses

Before testing our hypotheses with only religious people, we will filter only those that indicated that they are religious

# Filter data to include only religious participants (excludes "None", "Anti-religious", "Athiest")
data_rel <- data %>%
  filter(religious_affiliation %in% c("Christian", "Muslim", "Buddhist", 
                                      "Hindu", "Pagan", "Spiritual", "Spiritualist", "Deist"))

# Reversing the scoring for the SOI_a (saving it as SOI_a_r) and standardizing all variables
data_rel <- data_rel %>% 
  mutate(
    SOI_a_r = (max(data_rel$SOI_a) + 1) - SOI_a,
    SOI_a_r_z = scale(SOI_a_r),
    CRS_z = scale(CRS),
    religious_importance_z = scale(religious_importance),
    GENE_p_z = scale(GENE_p),
    GENE_p_2_z = scale(GENE_p_2),
    CONV_f_z = scale(CONV_f),
    TDDS_p_z = scale(TDDS_p)
  )

# Ensure that the variables are plain numeric vectors
data_rel$TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)
data_rel$CONV_f_z <- as.numeric(data_rel$CONV_f_z)
data_rel$CRS_z <- as.numeric(data_rel$CRS_z)
data_rel$SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
data_rel$religious_importance_z <- as.numeric(data_rel$religious_importance_z)
data_rel$GENE_p_z <- as.numeric(data_rel$GENE_p_z)
data_rel$GENE_p_2_z <- as.numeric(data_rel$GENE_p_2_z)

# Because the boot function was throwing a fit about this earlier
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)

Now, there is n = 141 participants for analysis here. All are religious.

Testing the Baseline Model

Before moving on, we will first test the baseline model where we regress the dependent variable (CRS) on the independent variable (TDDS_p), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions.

# Fitting the model
baseline_model_rel <- lm(CRS_z ~ TDDS_p_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(baseline_model_rel, which = 1)

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(baseline_model_rel))

   

   # QQ-plot for normality of residuals
   qqnorm(residuals(baseline_model_rel))
   qqline(residuals(baseline_model_rel), col = "red")

   

   • The histogram of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.

   • Because our inferential criteria do not require parametric tests, we should not be in too much trouble with a non-normality of residuals. We were relying on a parametric t-test to assess whether pathogen disgust is related to religiosity in the first place, however. Because we have failed the assumption of normality of residuals, we will use a bootstrapped test to make our inference here as well.

3. Multicollinearity: Not applicable becuase only one predictor

Summarizing the Model

# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(baseline_model_rel, col = "black", lwd = 1)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_baseline_model_rel <- Boot(baseline_model_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_baseline_model_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE    bootMed
(Intercept) 2.3086e-16 -7.2771e-05 0.079804 -0.0022716
TDDS_p_z    3.1293e-01  1.9212e-03 0.086146  0.3212187

confint(boot_baseline_model_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1471255 0.1619620
TDDS_p_z     0.1143230 0.4582468

# Removing the CRS_z and TDDS_p_z vectors in the environment that the Boot package required outside of the data.frame
rm(CRS_z, TDDS_p_z)

• We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and religiosity.

• Based on the estimates, the standardized coefficient for pathogen disgust is around .31. This is slightly higher than for the whole sample, including non-religious people (with was .24).

  • In addition, this relationship is significant, as the BCIs do not contain zero (upper = .47, lower = .15).

This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.

Hypothesis 1: Religious Subgroup Analyses

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religiosity.

# Fitting the models
H1_M1_rel <- lm(CONV_f_z ~ TDDS_p_z, data = data_rel)
H1_M2_rel <- lm(CRS_z ~ CONV_f_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_M1_rel, which = 1) # For model 1

   

   plot(H1_M2_rel, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_M1_rel)) # For model 1

   

   hist(residuals(H1_M2_rel)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H1_M1_rel))
   qqline(residuals(H1_M1_rel), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H1_M2_rel))
   qqline(residuals(H1_M2_rel), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just below 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$CONV_f_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Traditionalism",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M1_rel, col = "black", lwd = 1)

# Summarizing the model
summary(H1_M1_rel)

Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data_rel)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4199 -0.5498 -0.0864  0.6819  2.9068 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.454e-16  8.428e-02   0.000    1.000
TDDS_p_z    -7.504e-02  8.458e-02  -0.887    0.377

Residual standard error: 1.001 on 139 degrees of freedom
Multiple R-squared:  0.005631,  Adjusted R-squared:  -0.001523 
F-statistic: 0.7871 on 1 and 139 DF,  p-value: 0.3765

• We can see from the Scatter Plot that there is a slight negative relationship between pathogen disgust and traditionalism, but this relationship is not significant (B = -.07, t(139) = -.887, p = .377).

  • This negative result (while there is a null relationship in the population) could be an artifact of collider bias, where to be religious you could either be traditional or disgust sensitive and therefore in the population of religious people there is a slight negative relationship between these.

  • Although this precludes the possibility that adherence to tradition mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data_rel$CONV_f_z, data_rel$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Traditionalism",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H1_M2_rel, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data_rel$CONV_f_z)
CRS_z <- as.numeric(data_rel$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_rel <- Boot(H1_M2_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_M2_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE    bootMed
(Intercept) 2.6806e-16 -0.00086015 0.083509 -0.0012939
CONV_f_z    8.6179e-02  0.00724783 0.089536  0.0925694

confint(boot_H1_M2_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1622856 0.1598681
CONV_f_z    -0.1106986 0.2442643

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z)

• We can see from the Scatter Plot that there is a slight positive relationship between traditionalism and religiosity.

• Based on the estimates, the standardized coefficient for traditionalism is around .09.

  • However, this relationship is not significant, as the BCIs contain zero (upper = .25, lower = -.10).

Full Mediation Model

The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H1_full_rel <- lm(CRS_z ~ TDDS_p_z + CONV_f_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H1_full_rel, which = 1)

   

   • There is a trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H1_full_rel))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H1_full_rel))
   qqline(residuals(H1_full_rel), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H1_full_rel)

   TDDS_p_z CONV_f_z 
   1.005663 1.005663 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data_rel$CONV_f_z)
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_rel <- Boot(H1_full_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H1_full_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original    bootBias   bootSE   bootMed
(Intercept) 2.4689e-16 -9.7424e-05 0.079554 -0.001796
TDDS_p_z    3.2121e-01  4.5498e-04 0.088147  0.325711
CONV_f_z    1.1028e-01  6.8818e-03 0.078386  0.119082

confint(boot_H1_full_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.14873819 0.1608342
TDDS_p_z     0.11979446 0.4743443
CONV_f_z    -0.06180088 0.2456673

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results_rel <- mediate(model.m = H1_M1_rel,
        model.y = H1_full_rel,
        treat = "TDDS_p_z",
        mediator = "CONV_f_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(H1_mediate_results_rel)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value   
ACME           -0.00828     -0.03816         0.02   0.550   
ADE             0.32121      0.15388         0.48   0.002 **
Total Effect    0.31293      0.14954         0.47   0.002 **
Prop. Mediated -0.02644     -0.15190         0.08   0.548   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 141 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs not containing zero, pathogen disgust, but not traditionalism, is significantly positively related to religiosity.

  • The coefficients: for pathogen disgust (.32) and traditionalism (.11).

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = -.008, BCI = [-.044, .02]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, including when looking at religious people only.

Hypothesis 2: Religious Subgroup Analyses

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religiosity.

# Fitting the models
H2_M1_rel <- lm(GENE_p_z ~ TDDS_p_z, data = data_rel)
H2_M2_rel <- lm(CRS_z ~ GENE_p_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_M1_rel, which = 1) # For model 1

   

   plot(H2_M2_rel, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_M1_rel)) # For model 1

   

   hist(residuals(H2_M2_rel)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H2_M1_rel))
   qqline(residuals(H2_M1_rel), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H2_M2_rel))
   qqline(residuals(H2_M2_rel), col = "red")

   

   • Model 1:

     • The histogram and QQ-plot suggest that the residuals are relatively normally distributed for Model 1.

   • Model 2:

     • The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just below th mean and just below 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$GENE_p_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Out-Group Avoidance (GENE_p)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M1_rel, col = "black", lwd = 1)

# Summarizing the model
summary(H2_M1_rel)

Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data_rel)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0571 -0.6126 -0.0489  0.6251  3.8059 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.021e-16  8.438e-02   0.000    1.000
TDDS_p_z    5.742e-02  8.468e-02   0.678    0.499

Residual standard error: 1.002 on 139 degrees of freedom
Multiple R-squared:  0.003297,  Adjusted R-squared:  -0.003873 
F-statistic: 0.4598 on 1 and 139 DF,  p-value: 0.4988

• We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .057, t(139) = .678, p = .499).

  • Although this precludes the possibility that out-group avoidance mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.

Model 2

# Create a scatterplot with the regression line to visualize
plot(data_rel$GENE_p_z, data_rel$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Out-Group Avoidance (GENE_p)",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H2_M2_rel, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data_rel$GENE_p_z)
CRS_z <- as.numeric(data_rel$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_rel <- Boot(H2_M2_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_M2_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE     bootMed
(Intercept)  2.6622e-16 -0.00024138 0.084047  0.00042404
GENE_p_z    -5.4154e-02 -0.00716174 0.088882 -0.05820085

confint(boot_H2_M2_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1666378 0.1610149
GENE_p_z    -0.2223690 0.1094601

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z)

• We can see from the Scatter Plot that there is a negative trend in the relationship between out-group avoidance and religiosity.

  • Could be due to collider bias.

• Based on the estimates, the standardized coefficient for out-group avoidance is around -.05.

  • However, this relationship is not significant, as the BCI contains zero (upper = .12, lower = -.24).

Full Mediation Model

The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.

# Fitting the full mediation model
H2_full_rel <- lm(CRS_z ~ TDDS_p_z + GENE_p_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H2_full_rel, which = 1)

   

   • There is a slight trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H2_full_rel))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H2_full_rel))
   qqline(residuals(H2_full_rel), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H2_full_rel)

   TDDS_p_z GENE_p_z 
   1.003308 1.003308 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data_rel$GENE_p_z)
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_rel <- Boot(H2_full_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H2_full_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original    bootBias   bootSE    bootMed
(Intercept)  2.4549e-16 -7.8376e-05 0.080112 -0.0015804
TDDS_p_z     3.1709e-01  1.0441e-04 0.087176  0.3234525
GENE_p_z    -7.2362e-02 -4.9109e-03 0.079234 -0.0740654

confint(boot_H2_full_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %     97.5 %
(Intercept) -0.1502413 0.15947547
TDDS_p_z     0.1276348 0.46989918
GENE_p_z    -0.2260344 0.07345632

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results_rel <- mediate(model.m = H2_M1_rel,
        model.y = H2_full_rel,
        treat = "TDDS_p_z",
        mediator = "GENE_p_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H2_mediate_results_rel)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value   
ACME           -0.00416     -0.02253         0.02   0.708   
ADE             0.31709      0.14960         0.47   0.002 **
Total Effect    0.31293      0.14954         0.47   0.002 **
Prop. Mediated -0.01328     -0.08283         0.06   0.710   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 141 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.

  • The coefficients: for pathogen disgust (.32) and out-group avoidance (-.07).

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = -.004, BCI = [-.023, .02]).

  • This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even if only assessing religious participants.

Hypothesis 3: Religious Subgroup Analyses

Preliminary Simple Linear Regression Models

As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religiosity.

# Fitting the models
H3_M1_rel <- lm(SOI_a_r_z ~ TDDS_p_z, data = data_rel)
H3_M2_rel <- lm(CRS_z ~ SOI_a_r_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_M1_rel, which = 1) # For model 1

   

   plot(H3_M2_rel, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the either model.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_M1_rel)) # For model 1

   

   hist(residuals(H3_M2_rel)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(H3_M1_rel))
   qqline(residuals(H3_M1_rel), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(H3_M2_rel))
   qqline(residuals(H3_M2_rel), col = "red")

   

   • Both Models:

     • The QQ-plots for both models do not indicate that the residuals are normally distributed.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$SOI_a_r_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Pathogen Disgust",
     ylab = "Restricted Sociosexual Attitudes",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M1_rel, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_rel <- Boot(H3_M1_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M1_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original   bootBias   bootSE   bootMed
(Intercept) 1.3407e-15 0.00077784 0.081682 0.0028432
TDDS_p_z    1.9972e-01 0.00086515 0.080392 0.1985761

confint(boot_H3_M1_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1664251 0.1462323
TDDS_p_z     0.0464425 0.3473499

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)

• We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.

• Based on the estimates, the standardized coefficient for pathogen disgust is around .19.

  • In addition, this relationship is significant, as the BCI excludes zero (upper = .35, lower = .03).

Model 2

# Create a scatterplot with the regression line to visualize
plot(data_rel$SOI_a_r_z, data_rel$CRS_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Restricted Sociosexual Attitudes",
     ylab = "Religiosity (CRS)",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(H3_M2_rel, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
CRS_z <- as.numeric(data_rel$CRS_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_rel <- Boot(H3_M2_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_M2_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE    bootMed
(Intercept) -2.6447e-16 -0.0013189 0.078758 -0.0016614
SOI_a_r_z    3.8309e-01 -0.0043342 0.080185  0.3809897

confint(boot_H3_M2_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                 2.5 %    97.5 %
(Intercept) -0.1610559 0.1580375
SOI_a_r_z    0.2218491 0.5358196

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z)

• We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.

• Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .38.

  • In addition, this relationahip is significant, as the BCI excludes zero (upper = .54, lower = .22).

Full Mediation Model

The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.

# Fitting the full mediation model
H3_full_rel <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z, data = data_rel)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(H3_full_rel, which = 1)

   

   • There is a slight trend towards more negative residual values at the higher fitted values of the model. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(H3_full_rel))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(H3_full_rel))
   qqline(residuals(H3_full_rel), col = "red")

   

   • The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just below the mean and just below 1 SD.

     • Because of this failure of this assumption, we will conduct a bootstrapped test of the predictors for the model.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(H3_full_rel)

    TDDS_p_z SOI_a_r_z 
    1.041546  1.041546 

   • The VIF is good here.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_rel <- Boot(H3_full_rel, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_H3_full_rel) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE    bootMed
(Intercept) -2.1680e-16 -0.0012018 0.076513 -0.0015274
TDDS_p_z     2.4625e-01  0.0023376 0.086831  0.2520226
SOI_a_r_z    3.3391e-01 -0.0043978 0.082107  0.3283799

confint(boot_H3_full_rel) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.14622859 0.1585309
TDDS_p_z     0.05404243 0.4148393
SOI_a_r_z    0.17970551 0.4927243

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results_rel <- mediate(model.m = H3_M1_rel,
        model.y = H3_full_rel,
        treat = "TDDS_p_z",
        mediator = "SOI_a_r_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results of the mediation analysis
summary(H3_mediate_results_rel)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value   
ACME             0.0667       0.0141         0.14   0.012 * 
ADE              0.2462       0.0806         0.41   0.006 **
Total Effect     0.3129       0.1495         0.47   0.002 **
Prop. Mediated   0.2131       0.0479         0.53   0.014 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 141 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs, both pathogen disgust and restricted sociosexual attitudes are significantly positively related to religiosity, with the unique effect of restricted sociosexual attitudes being substantially larger (B = .334, BCI =[.173, .501]) than the unique effect of pathogen disgust (B = .246, BCI = [.083, .403]).

• Mediation Analysis Results:

  • The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .067, BCI = [.013, .13]).

    • Given this result, we have evidence for the weak version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is partially mediated by a monogamous mating strategy, including when only religious individuals are included in the analysis.

  • The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .21. That is, 21% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.044, .54]).

    • Given this result, we do not have evidence for the strong version of Hypothesis 3. That is, the relationship between pathogen disgust and religiosity is not fully mediated by a monogamous mating strategy, including when only religious individuals are included.

Ultimately, the results of including only religious individuals in the analysis support the results thus far. They are in line with the weak version of the sexual strategies hypothesis only. In addition, the proportion mediated is slightly smaller when analyzing only religious individuals (21% vs. 34%). This suggests that while the sexual strategies hypothesis may explain some differences in religiosity between religious people, it is better at explaining individual differences in religiosity including those that are not explicitly religious. This may simply be a feature of restricted range, or it may suggest that the sexual strategies hypothesis can explain both why religious people differ in their “devotion” and why some people are religious and others are not.

Sexual Disgust, Sociosexual Attitudes, and Traditionalism (SST Model)

BIS reactivity may encourage traditional cultural values (like religion) as a means of mitigating social behaviors (like sex) that are more likely to result in the spread of infectious disease. We have focused thus far on religion in this analysis as the cultural value influenced by BIS reactivity and the specific mechanisms through which this may be the case. We will now analyze whether sociosexual attitudes mediate the relation between sexual disgust and traditionalism.

Testing the Baseline Model

We will first test the baseline model where we regress the dependent variable (CONV_f_z) on the independent variable (TDDS_s_z), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions. We will first need to create the TDDS_s_z variable by standardizing TDDS_s.

# Creating the standardized sexual disgust variable
data$TDDS_s_z <- scale(data$TDDS_s)
  # Ensure it is numeric so the boot function will run
data$TDDS_s_z <- as.numeric(data$TDDS_s_z)

# Fitting the model
baseline_model_SST <- lm(CONV_f_z ~ TDDS_s_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(baseline_model_SST, which = 1)

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(baseline_model_SST))

   

   # QQ-plot for normality of residuals
   qqnorm(residuals(baseline_model_SST))
   qqline(residuals(baseline_model_SST), col = "red")

   

   • The histogram of the residuals of the model indicate that the residuals are normally distributed.

3. Multicollinearity: Not applicable becuase only one predictor

Summarizing the Model

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_s_z, data$CONV_f_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Sexual Disgust",
     ylab = "Traditionalism",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(baseline_model_SST, col = "black", lwd = 1)

# Summarizing the model
summary(baseline_model_SST)

Call:
lm(formula = CONV_f_z ~ TDDS_s_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3762 -0.5975  0.0059  0.6273  2.7846 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -6.880e-17  5.842e-02   0.000   1.0000  
TDDS_s_z     1.307e-01  5.852e-02   2.234   0.0263 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9931 on 287 degrees of freedom
Multiple R-squared:  0.01709,   Adjusted R-squared:  0.01366 
F-statistic: 4.989 on 1 and 287 DF,  p-value: 0.02628

• We can see from the Scatter Plot that there is a slight positive relationship between sexual disgust and religiosity.

• Based on the estimates, the (significant) standardized coefficient for sexual disgust is around .131 (p = .026).

This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.

Preliminary Simple Linear Regression Models

As mentioned above, to test our mediaiton model, we must first establish that sexual disgust affects sociosexual attitudes (i.e., SOI_a_r_z) and that sociosexual attitudes influences traditionalism.

# Fitting the models
SST_M1 <- lm(SOI_a_r_z ~ TDDS_s_z, data = data)
SST_M2 <- lm(CONV_f_z ~ SOI_a_r_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(SST_M1, which = 1) # For model 1

   

   plot(SST_M2, which = 1) # For model 2

   

   • The residuals vs. fitted plot does not stand out as violating the assumption of linearity or heteroscedasticity for the second model, but there is a trend towards more negative residual values at the higher fitted values of the model for model 1. This may be problematic.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(SST_M1)) # For model 1

   

   hist(residuals(SST_M2)) # For model 2

   

   # QQ-plot for normality of residuals of Model 1
   qqnorm(residuals(SST_M1))
   qqline(residuals(SST_M1), col = "red")

   

   # QQ-plot for normality of residuals of Model 2
   qqnorm(residuals(SST_M2))
   qqline(residuals(SST_M2), col = "red")

   

   • Both Models:

     • The histogram and QQ-plot suggest that the residuals are relatively normally distributed for both models.

   • We will conduct bootstrapping for Model 1 because of the ambivalent residuals plot above.

3. Multicollinearity: Not applicable because each model only has one predictor.

Summarizing the Models

Model 1

# Create a scatterplot with the regression line to visualize
plot(data$TDDS_s_z, data$SOI_a_r_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Restricted Sociosexual Attitudes",
     ylab = "Sexual Disgust",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(SST_M1, col = "black", lwd = 1)

# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_s_z <- as.numeric(data$TDDS_s_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_SST_M1 <- Boot(SST_M1, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_SST_M1) # Getting the coefficients

Number of bootstrap replications R = 1000 
               original   bootBias   bootSE   bootMed
(Intercept) -1.5856e-15 -0.0020961 0.046731 -0.001715
TDDS_s_z     5.8840e-01  0.0002019 0.042232  0.588610

confint(boot_SST_M1) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %     97.5 %
(Intercept) -0.09047559 0.09164489
TDDS_s_z     0.50214194 0.66703586

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_s_z)

• We can see from the Scatter Plot that there is a strong positive relationship between sexual disgust and restricted sociosexual attitudes (B = .588, BCI = [.491, .679]).

Model 2

# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$CONV_f_z,
     main = "Scatter Plot with Linear Regression Line",
     xlab = "Restricted Sociosexual Attitudes",
     ylab = "Traditionalism",
     pch = 16, # For solid circle points
     col = "black")
# Add regression line
abline(SST_M2, col = "black", lwd = 1)

# Summarizing the model
summary(SST_M2)

Call:
lm(formula = CONV_f_z ~ SOI_a_r_z, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5045 -0.6536  0.0191  0.6226  2.7662 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.739e-16  5.826e-02   0.000   1.0000  
SOI_a_r_z   1.500e-01  5.836e-02   2.571   0.0107 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9904 on 287 degrees of freedom
Multiple R-squared:  0.02251,   Adjusted R-squared:  0.0191 
F-statistic: 6.609 on 1 and 287 DF,  p-value: 0.01065

• We can see from the Scatter Plot that there is a slight positive relationship between sociosexual attitudes and traditionalism.

  • This relationship is significant (B = .15, t(287) = 2.571, p = .011).

Full Mediation Model

# Fitting the full mediation model
SST_full <- lm(CONV_f_z ~ TDDS_s_z + SOI_a_r_z, data = data)

Assumptions

1. Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values

   # Plotting the residual vs. fitted values
   plot(SST_full, which = 1)

   

   • The plot looks good.

2. Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)

   # Histogram for normality of residuals
   hist(residuals(SST_full))

   

   # QQ-plot for normality of residuals of the full model
   qqnorm(residuals(SST_full))
   qqline(residuals(SST_full), col = "red")

   

   • The histogram looks good.

3. Multicollinearity: Variance Inflation Factor (VIF) from car package

   # Calculating VIF
   vif(SST_full)

    TDDS_s_z SOI_a_r_z 
    1.529543  1.529543 

   • The VIF is higher in this model than all previous models due to the correlation between sexual disgust and restricted sociosexual attitudes. However, it is within an acceptable range for our analysis.

Summarizing and Testing for Partial/Full Mediation

# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
TDDS_s_z <- as.numeric(data$TDDS_s_z)
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)

# Setting the seed
set.seed(1042557)

# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_SST_full <- Boot(SST_full, R = 1000) # 1,000 resamples

# Summarizing the estimates and CIs
summary(boot_SST_full) # Getting the coefficients

Number of bootstrap replications R = 1000 
              original   bootBias   bootSE     bootMed
(Intercept) 1.0852e-16 -0.0016315 0.059488 -0.00096288
TDDS_s_z    6.4909e-02 -0.0003628 0.071978  0.06578146
SOI_a_r_z   1.1184e-01  0.0030682 0.072723  0.11451771

confint(boot_SST_full) # Getting the BCIs

Bootstrap bca confidence intervals

                  2.5 %    97.5 %
(Intercept) -0.11907106 0.1159598
TDDS_s_z    -0.07670702 0.1985219
SOI_a_r_z   -0.03639622 0.2529179

# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, TDDS_s_z, SOI_a_r_z)

# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
SST_mediate_results <- mediate(model.m = SST_M1,
        model.y = SST_full,
        treat = "TDDS_s_z",
        mediator = "SOI_a_r_z",
        boot = TRUE, # Bootstrap the CI for the indirect effect
        sims = 1000) # Do 1000 resamples

Running nonparametric bootstrap

# Printing the results for the mediation analysis
summary(SST_mediate_results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value  
ACME             0.0658      -0.0168         0.14    0.13  
ADE              0.0649      -0.0626         0.20    0.39  
Total Effect     0.1307       0.0218         0.25    0.03 *
Prop. Mediated   0.5034      -0.3460         2.22    0.16  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 289 


Simulations: 1000 

• Multiple Regression Results:

  • Based on the bootstrapped CIs containing zero, neigher sexual disgust nor restricted sociosexual attitudes are significantly positively related to traditionalism.

  • The coefficients for sexual disgust (.065) and sociosexual attitudes (.112) are quite small.

• Mediation Analysis Results:

  • The ACME indicates that there is no significant indirect effect of sexual disgust on traditionalism through restricted sociosexual atitudes (B = .066, BCI = [-.02, .15]).

Even though the proportion mediated is high (50%), it seems that the lack of a relationship between the predictors (over and above one another) and traditionalism means that this large estimated proportion mediated is not reliable.

Relating Pathogen Disgust and Religiosity After Removing Regional Differences

To attempt to remove the effect of national levels of religiosity and disgust sensitivity from influencing the relationship between disgust and religiosity at the individual level, we will see if we can isolate the individual effect by running a multi-level model with random intercepts. This will account for mean levels of religiosity and isolate individual-level effects from their confounding influence. However, we do not have enough participants in most countries to do this at the country level. Therefore, I will split the nationalities into world region groups. This is admittedly a compromise move, because using national-level groupings would be better. In addition, for many nations, there is not an obvious world region given a lack of representation of participants in countries of those regions. For example, I will exclude participants in areas around the Asia Pacific and participants in South America. We have such nationalities in this sample, but there are too few participants in these regions to justify including them. The grouping that I could come up with that maximizes the number of people in each group while retaining regional coherency was a three-group division: (1) Europe (n = 104), (2) North America (n = 86), and (3) Africa and the Middle East (n = 59). We will split participants up into these regions based on their nationality and fit a multi-level model with a random intercept of region. This will allow us to isolate the individual influence of disgust and religiosity without these three regions’ influence on religiosity scores.

# Load necessary packages
library(lme4) # For fitting the multi-level model
library(sjPlot)      # For visualization

Learn more about sjPlot with 'browseVignettes("sjPlot")'.

# Assign regions based on nationality
MLM_data_region <- data %>%
  mutate(region = case_when(
    nationality %in% c("Canada", "United States", "Mexico") ~ "North America",
    nationality %in% c("Belgium", "Austria", "Czech Republic", "Finland", "France", 
                       "Germany", "Greece", "Hungary", "Ireland", "Italy", "Poland", 
                       "Portugal", "Romania", "Slovakia", "Spain", "Sweden", 
                       "United Kingdom") ~ "Europe",
    nationality %in% c("Algeria", "Egypt", "Eritrea", "Ghana", "Kenya", 
                       "South Africa", "Saudi Arabia", "Tunisia", "Zimbabwe") ~ "Africa & Middle East",
    TRUE ~ NA_character_
  ))

MLM_data_region <- MLM_data_region %>% 
  filter(!is.na(region))

# Fit the multilevel model with standardized variables and region as random intercept
MLM_disg_rel_region <- lmer(
  CRS_z ~ TDDS_p_z + (1 | region),
  data = MLM_data_region,
  REML = FALSE
)

# Fixed effects plot
sjPlot::plot_model(MLM_disg_rel_region, type = "pred", terms = "TDDS_p_z")

# Random effects plot
sjPlot::plot_model(MLM_disg_rel_region, type = "re")

# Faceted scatter plot
ggplot(MLM_data_region, aes(x = TDDS_p_z, y = CRS_z)) +
  geom_point(alpha = 0.3, color = "darkblue") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ region) +
  theme_minimal() +
  labs(title = "Standardized Disgust Sensitivity vs. Religiosity Across Regions",
       x = "Standardized Disgust Sensitivity",
       y = "Standardized Religiosity")

`geom_smooth()` using formula = 'y ~ x'

• Looking at the fixed effects plot, we can see that the general tendency is for disgust sensitivity to positively predict religiosity.

• Looking at the random intercepts, we can see that while North America and Europe are similar in their levels of religiosity at their mean of disgust sensitivity, Africa and the Middle East is quite high. It is about .75 SD above the mean whereas Europe and North America are about .7 and .3 SD below the mean.

  • This may represent the concern of Reviewer 2, where countries with higher mean levels of religiosity and disgust sensitivity (e.g. those in the Middle East and Africa) may bias the individual-level results.

• Looking at the faceted scatter plots, it looks like there is a far more distinct positive relationship in Europe—and to a lesser degree in North America—than in Africa and the Middle East between disgust sensitivity and religiosity.

  • This may partially be due to restricted range, as the relationship within Africa and the Middle East may be constricted by not having enough range in the religiosity scale towards the higher end. The points seem to cluster towards the top of the graph for Africa and the Middle East.

Before moving on to summarizing and testing our model parameter, we will first test some assumptions.

Assumptions

Here we will test the assumptions of homogeneity of variance and normality of residuals for the model by plotting the residuals vs. the fitted values and a histogram and QQ-plot.

# Plotting the residual vs. fitted values
plot(MLM_disg_rel_region, which = 1)

# Histogram for normality of residuals
hist(residuals(MLM_disg_rel_region))

# QQ-plot for normality of residuals
qqnorm(residuals(MLM_disg_rel_region))
qqline(residuals(MLM_disg_rel_region), col = "red")

• The residuals vs. fitted values shows two clusters of residuals at the high and low and of the scale, and they are consistent with large differences in religiosity in our three regions. That is, the cluster on the high end of the scale should be Africa and the Middle East. The cluster on the low end of the scale should represent the United States and Europe. This may be a problem for inference, so we will conduct a bootstrapped test of the fixed effect of disgust sensitivity on religiosity.

• The histogram and QQ-plot look pretty good.

Model Inference

Now we will conduct a bootstrapped test of the fixed effect of disgust sensitivity on religiosity.

# Set seed for reproducibility
set.seed(1042557)

# Function to extract fixed effects
fixef_fun <- function(model) {
  fixef(model)
}

# Perform bootstrapping
boot_results_MLM_disg_rel_region <- bootMer(
  MLM_disg_rel_region,
  FUN = fixef_fun,
  nsim = 1000,
  type = "parametric",
  cl = cl,
  seed = 1042557
)

# Calculate 95% confidence intervals
boot_conf_interval_MLM_disg_rel_region <- apply(boot_results_MLM_disg_rel_region$t, 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)

# Convert to a more readable format
boot_conf_interval_MLM_disg_rel_region <- as.data.frame(boot_conf_interval_MLM_disg_rel_region)
rownames(boot_conf_interval_MLM_disg_rel_region) <- c("2.5%", "97.5%")

# View the estimates and BCI for fixed effect
summary(MLM_disg_rel_region)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: CRS_z ~ TDDS_p_z + (1 | region)
   Data: MLM_data_region

     AIC      BIC   logLik deviance df.resid 
   608.1    622.1   -300.0    600.1      245 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6269 -0.6843 -0.1220  0.6094  2.7162 

Random effects:
 Groups   Name        Variance Std.Dev.
 region   (Intercept) 0.4390   0.6626  
 Residual             0.6207   0.7878  
Number of obs: 249, groups:  region, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.21242    0.38598   0.550
TDDS_p_z     0.09887    0.05160   1.916

Correlation of Fixed Effects:
         (Intr)
TDDS_p_z -0.011

print(boot_conf_interval_MLM_disg_rel_region)

      (Intercept)     TDDS_p_z
2.5%   -0.5445795 -0.000325435
97.5%   0.9450968  0.199948070

• We can see that the random intercept of region accounts for a large amount of variance in religiosity. That is, 44% of the variance in religiosity is attributable to the mean level of religiosity within the region.

• However, we can also see that individual disgust sensitivity still is a very marginally significant predictor of religiosity, albeit a smaller size (B = .100, BCI = [.200, -.00003].

Ultimately this analysis suggests that while the relationship between disgust and religiosity may be non-zero, it is being inflated by the presence of Africa and the Middle East in the sample. Similarly, however, the relationship between disgust sensitivy and religiosity appears to be closer to B = 0 within Africa and the Middle East due to restricted range. This suggests that the relationship in Europe and North America may be higher even if accounting for regional variation. To explore the data further, I will fit a similar model, but this time only with North America and Europe as regions.

Relating Pathogen Disgust and Religiosity in the West

# Assign regions based on nationality
MLM_data_west <- data %>%
  mutate(region = case_when(
    nationality %in% c("Canada", "United States", "Mexico") ~ "North America",
    nationality %in% c("Belgium", "Austria", "Czech Republic", "Finland", "France", 
                       "Germany", "Greece", "Hungary", "Ireland", "Italy", "Poland", 
                       "Portugal", "Romania", "Slovakia", "Spain", "Sweden", 
                       "United Kingdom") ~ "Europe",
    TRUE ~ NA_character_
  ))

MLM_data_west <- MLM_data_west %>% 
  filter(!is.na(region))

# Fit the multilevel model with standardized variables and region as random intercept
MLM_disg_rel_west <- lmer(
  CRS_z ~ TDDS_p_z + (1 | region),
  data = MLM_data_west,
  REML = FALSE
)

# Fixed effects plot
sjPlot::plot_model(MLM_disg_rel_west, type = "pred", terms = "TDDS_p_z")

# Random effects plot
sjPlot::plot_model(MLM_disg_rel_west, type = "re")

# Faceted scatter plot
ggplot(MLM_data_west, aes(x = TDDS_p_z, y = CRS_z)) +
  geom_point(alpha = 0.3, color = "darkblue") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ region) +
  theme_minimal() +
  labs(title = "Standardized Disgust Sensitivity vs. Religiosity Across Regions",
       x = "Standardized Disgust Sensitivity",
       y = "Standardized Religiosity")

`geom_smooth()` using formula = 'y ~ x'

• Looking at the fixed effects plot, we can see that the general tendency is for disgust sensitivity to positively predict religiosity.

• Looking at the random intercepts, we can see that while North America and Europe are similar in their levels of religiosity at their mean of disgust sensitivity, with North America being slightly more religious.

• Looking at the faceted scatter plots,the relationships look similar to the previous plots. There appears to be a positive relationship in both Europe and North America.

Before moving on to summarizing and testing our model parameter, we will first test some assumptions.

Assumptions

Here we will test the assumptions of homogeneity of variance and normality of residuals for the model by plotting the residuals vs. the fitted values and a histogram and QQ-plot.

# Plotting the residual vs. fitted values
plot(MLM_disg_rel_west, which = 1)

# Histogram for normality of residuals
hist(residuals(MLM_disg_rel_west))

# QQ-plot for normality of residuals
qqnorm(residuals(MLM_disg_rel_west))
qqline(residuals(MLM_disg_rel_west), col = "red")

• The residuals vs. fitted values doesn’t show much of an issue. There may be more extreme positive residual values and more negative residual values in general. There may also be a slight tendency for residuals to be more negatively extreme towards the positive end of the fitted values. However, there are no glaring patterns indicating a failure of assumptions.

• The histogram and QQ-plot look pretty skewed, so we will do a bootstrap inference for the parameter.

Model Inference

Now we will conduct a bootstrapped test of the fixed effect of disgust sensitivity on religiosity.

# Set seed for reproducibility
set.seed(1042557)

# Function to extract fixed effects
fixef_fun <- function(model) {
  fixef(model)
}

# Perform bootstrapping
boot_results_MLM_disg_rel_west <- bootMer(
  MLM_disg_rel_west,
  FUN = fixef_fun,
  nsim = 1000,
  type = "parametric",
  cl = cl,
  seed = 1042557
)

# Calculate 95% confidence intervals
boot_conf_interval_MLM_disg_rel_west <- apply(boot_results_MLM_disg_rel_west$t, 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)

# Convert to a more readable format
boot_conf_interval_MLM_disg_rel_west <- as.data.frame(boot_conf_interval_MLM_disg_rel_west)
rownames(boot_conf_interval_MLM_disg_rel_west) <- c("2.5%", "97.5%")

# View the estimates and BCI for fixed effect
summary(MLM_disg_rel_west)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: CRS_z ~ TDDS_p_z + (1 | region)
   Data: MLM_data_west

     AIC      BIC   logLik deviance df.resid 
   468.4    481.4   -230.2    460.4      186 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4893 -0.6903 -0.2154  0.4676  2.5404 

Random effects:
 Groups   Name        Variance Std.Dev.
 region   (Intercept) 0.008748 0.09353 
 Residual             0.654947 0.80929 
Number of obs: 190, groups:  region, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.25509    0.08890  -2.869
TDDS_p_z     0.12206    0.06018   2.028

Correlation of Fixed Effects:
         (Intr)
TDDS_p_z 0.090 

print(boot_conf_interval_MLM_disg_rel_west)

      (Intercept)    TDDS_p_z
2.5%  -0.44074157 0.008464031
97.5% -0.07839034 0.240534554

• We can see that the random intercept of region accounts for very little variance in religiosity in this model. Only 0.8% of the variance in religiosity is attributable to the mean level of religiosity within the region.

• We can also see that individual disgust sensitivity now a slightly larger predictor of religiosity and is significant (B = .122, BCI = [.241, .008].

This analysis suggests that the relationship between disgust and religiosity in Europe and North America is more robust. This may suggest a relationship between disgust and religiosity that is specific to Christian cultures, because Europe and North America tend to be Christian.

Saving Data

Now we will save the final version of the data frame (data) as a data/final_data.csv.

# Saving the data
write.csv(data, "./data/final_data.csv")

References

Baron, R. M., & Kenny, D. A. (1986). The moderator–mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations. Journal of Personality and Social Psychology, 51(6), 1173–1182. https://doi.org/10.1037/0022-3514.51.6.1173

Beribisky, N., Mara, C. A., & Cribbie, R. A. (2020). An equivalence testing approach for evaluating substantial mediation. The Quantitative Methods for Psychology, 16(5), 424–441. https://doi.org/10.20982/tqmp.16.4.p424

Duckitt, J., Bizumic, B., Krauss, S. W., & Heled, E. (2010). A Tripartite Approach to Right-Wing Authoritarianism: The Authoritarianism-Conservatism-Traditionalism Model. Political Psychology, 31(5), 685–715. https://doi.org/10.1111/j.1467-9221.2010.00781.x

Dunwoody, P. T., & Funke, F. (2016). The Aggression-Submission-Conventionalism Scale: Testing a New Three Factor Measure of Authoritarianism. Journal of Social and Political Psychology, 4(2), Article 2. https://doi.org/10.5964/jspp.v4i2.168

Huber, S., & Huber, O. W. (2012). The Centrality of Religiosity Scale (CRS). Religions, 3(3), Article 3. https://doi.org/10.3390/rel3030710

Neuliep, J. W. (2002). Assessing the Reliability and Validity of the Generalized Ethnocentrism Scale. Journal of Intercultural Communication Research, 31(4), 201. https://twu.idm.oclc.org/login?url=https://search.ebscohost.com/login.aspx?direct=true&db=ufh&AN=10966365&site=eds-live&scope=site

Neuliep, J. W., & McCroskey, J. C. (1997). The development of a US and Generalized Ethnocentrism Scale. Communication Research Reports, 14(4), 385–398. https://doi.org/10.1080/08824099709388682

Olatunji, B., Adams, T., Ciesielski, B., David, B., Sarawgi, S., & Broman-Fulks, J. (2012). The Three Domains of Disgust Scale. Assessment, 19, 205–225. https://doi.org/10.1177/1073191111432881

Penke, L., & Asendorpf, J. B. (2008). Beyond global sociosexual orientations: A more differentiated look at sociosexuality and its effects on courtship and romantic relationships: Journal of Personality and Social Psychology. Journal of Personality and Social Psychology, 95(5), 1113–1135. https://doi.org/10.1037/0022-3514.95.5.1113

Simpson, J. A., & Gangestad, S. W. (1991). Individual differences in sociosexuality: evidence for convergent and discriminant validity. Journal of Personality and Social Psychology, 60(6), 870–883. https://doi.org/10.1037//0022-3514.60.6.870

Tybur, J. M., Lieberman, D. L., & Griskevicius, V. G. (2009). Microbes, mating, and morality: Individual differences in three functional domains of disgust. Journal of Personality and Social Psychology, 29, 103–122. https://doi.org/10.1037/a0015474

Yu, Z., Bali, P., Tsikandilakis, M., & Tong, E. M. W. (2022). “Look not at what is contrary to propriety”: A meta-analytic exploration of the association between religiosity and sensitivity to disgust. British Journal of Social Psychology, 61(1), 276–299. https://doi.org/10.1111/bjso.12479